On Fitting Combined Integral and Differential Reaction Kinetic Data. Rate Modeling for Catalytic Hydrogenation of Butadiene Pradeep P. Sane,l Roger E. Eckert, and John M. Woods* School of Chemical Engineering, Purdue University, Lafayette, lnd. 47907
A regression ’technique for the analysis of combined integral and differential kinetic data is developed by the correct use of average rate of reaction as a response variable. The technique is used to analyze data on the hydrogenation of butadiene in terms of each of two plausible Hougen-Watson rate models. A simple lack-of-fit test discriminates between the two models and establishes the adequacy of one.
Introduction
Experimental Results
In the mechanistic modeling of heterogeneous catalytic reaction rates the basic problems involved are: (1) identification of an adequate model and (2) estimation of the parameters within that model. Hougen and Watson (1943) have pioneered techniques to formulate mechanistic models for the gas-solid catalytic reactions. For the same chemical reaction different assumptions concerning various steps of the reaction lead to different models. Advanced statistical techniques of designing experiments (e.g., Box and Hill, 1967; Froment and Mezaki, 1970; Hill, et al., 1968; Mezaki and Happel, 1969) are powerful tools for model discrimination and precise parameter estimation. These techniques aim mainly a t minimizing the amount of experimentation. If extensive data over wide ranges of experimental conditions are available, even simple techniques like linear and nonlinear regression analysis sometimes can be used effectively. Complete reaction kinetic studies often include both integral and differential data. Statistical methods for handling either of these classes of data to obtain the bkst parameter estimates for a reaction rate model have been reviewed by Kittrell (1970). The correct use of an appropriate response variable is of vital importance in the analysis of the corn bined integral and differential data. Recently, Sane, et al. (1972), outlined a method to analyze the combined data. In the present study, we shall apply this method to the modeling of butadiene hydrogenation rates. This reaction has been studied by Gerheim (1966) and Kahney (1969) in a packed bed flow reactor using platinum on an alumina support as a catalyst. The hydrogenation of butadiene may be represented schematically by the equation
A total of 145 data points summarized in Table I are available for the analysis. As an example of the distribution of points in p~~ and P ~ Dspace the data a t approximately 367°K is shown in Figure 1. The fractional conversion of butadiene was very low ( E 1%) in the region investigated by Gerheim. Extremely low butadiene partial pressures studied by Kahney led to high fractional conversion. In Kahney’s region of investigation, the amount of diluent (helium gas) used was large compared to the amount of butadiene. Also changes in the total number of moles and the total pressure were insignificant; consequently, the approximation
H
n
LSbutadiene 4 butenes A butane
(1)
Chromatographic analysis of the feed and product gases was used to determine the fraction of the butadiene feed hydrogenated in the reactor. After some preliminary analyses of the different subsets of the data, Kahney has recommended two plausible Hougen-Watson models for these data. First we shall present briefly the available data on the hydrogenation of butadiene and the two proposed rate models. The formulation of regression models for the analysis of combined data follows. Finally, both of the proposed models are analyzed and a simple lack o f f i t test is performed to decide the superior model. l
Present address, Polyolefins Industries Ltd., P.O. Box 1151, Bombay
1 , India. 52
Ind. Eng. Chern., Fundarn., Vol. 13, No. 1, 1974
p’= nBD(” - nsJf) nBD(I)
-
pBD(l)
- pBdf’
PBD(’)
(2)
holds for the analysis to follow. The temperature of the bulk gas a t various points along the length of the reactor was measured during the experimentation. The temperature rise in the reactor was less than 1°K for most of the data points. Intraparticle temperature gradients were estimated and found to be negligible. Intraparticle concentration gradients and interparticle temperature gradients were studied by simulating the laboratory reactor (Sane, 1972) and did not change the results appreciably. Moreover, the residual error in the regression analysis was nearly the same when including or excluding such an effect. For example, 26 of Kahney’s data points a t 357°K showed the residual error to be 2.16 x g mole/(sec) (g of catalyst) when the intraparticle diffusion effect was neglected, as against 2.156 X when the effect was included. Therefore, these effects were ignored in the regression model. The activity of the catalyst varied during the period of experimentation. A standard activity was defined and relative activities a t various times were determined by frequently making “standard runs” during the investigation period. The details of this procedure are outlined by both Gerheim and Kahney. This procedure enabled them to assign an approximate value of the catalyst activity to each data point. The value of the rate constant was modified by the catalyst activity as (3) h’ = h x (activity) Kinetic Models
Initial fitting of power law models to the various subsets of the data and the literature study of the Hougen-Watson models postulated for the hydrogenation reactions led
Table I. Summary of Experimental Data No. of
Investigator
data points
P H ~ ,mm
Gerheim Kahney
97 48
176.6-513 .6 215 .6-388.9
500t QlI 400t@
PBD,
mm
33-264 2 .28-18.25
a@a@
00
IO0
c
.
THE NUMBERS IN THE CIRCLE REPRESENT THE NUMBER OF
EXPERIMENTS PERFORMED I N EACH S M A L L REGION
PBO
Figure 1. Experiments on 367°K
PED
mm H g
-
space at approximately
PH?
Kahney to propose the following two models as the plausible ones. Model I
R=
kpH21'5pBD
(1 + KBDPBD + KBEPBE)'
(4)
Model I1
PHE,
T , OK
mm
...
363.1-387.1 357.3, 367.3
0, 51, 102
Regression Model for Combined Data The major problem in fitting the combined integral and differential kinetics data is the choice of an appropriate response variable. Sane, e t al. (1972), have shown that the average rate of reaction, as defined by eq 7 in the present case, is such a variable. -
KBD,
k = k, exp(-A/T) KBD
= K B D , exp(-B/T)
R-
KBE
= KB,
exp(-C/T)
-
The rate models given by eq 4 and 5 are obtained by considering the following sequence of reaction steps.
+s
C4big)
+ 2s
H,'"
+
(C&)s (C,H,)ss
+
-
(C4H6)S
(a)
2(H),
(b)
Hs + (C4H7)s
(C4H&
(C&)s
C4&'gJ
+ (H),
+S
(C)
(4 (e)
Step d in the forward direction is assumed to be rate controlling. The reverse of step d is neglected and all the rest of the steps are considered to be in a quasiequilibrium state. Further, the adsorption coefficient of hydrogen is considered to be very small compared to those of butadiene and butenes. The single site adsorption of butadiene shown in step a leads to Model I, while the dual site adsorption shown below leads to Model 11. C4Hig'
+ 2s
==
(C,H&
(a')
The earlier workers fitted both the models to the differential data collected a t high partial pressures of butadiene
~(~)p~,"'&
22,400PW
(7)
Under differential conditions, the average rate and initial rate are nearly equivalent.
R'"
(6)
* 1%
up to 95%
(>30 mm). In this application data were considered to be suitable for differential anal$& if the fractional conversion of butadiene was about 1% or less. The results of nonlinear regression analyses by Kahney are shown in Table 11. The value of the coefficient of multiple determination, R2, a t 375°K is the same for both the models. Model I gives a higher R2 a t 365°K but the temperature dependency of k in Model I is abnormal since the k is seen t o decrease with the increase in temperature. Therefore, analysis of the differential data alone is insufficient to discriminate between Models I and 11. Data at the higher partial pressures of butadiene neglected a region critical to model discrimination. Only in the region of low partial pressure of butadiene did the difference between rates predicted by the two models become large. In this region of experimentation a small absolute conversion is a large fractional conversion of the initial butadiene. Differential analysis is valid only for small fractional conversions. Thus integral. analysis was necessary for the low partial pressure data. Finally, the integral data must be combined with the differential data to decide the superior rate model.
(5) The temperature dependency of the parameters k , and K Hcan ~ be taken into account by the following.
Level of butadiene conversion
M(T"', PH?''),PBD(I)jP B E ' I ) )
(8) M is the functional form of the rate model under consideration. For integral data the near equivalence of the average and the initial rates breaks down because of the relatively large ohanges in the amount present of a t least one of the species, butadiene in this case. Under these conditions the fractional conversion usually is taken to be the response variable. When differential and integral data are combined, the conversion is rarely a good response variable. The standard deviation of the conversion for integral data is often of the same magnitude as the conversion itself for the differential data; the standard deviation for the differential conversion is much lower. Furthermore, it is reasonable to assume that n B D ' l ' has a negligible vari) . according to eq 2 ance as compared to ~ B D ( ~Therefore the variance of x ' ~ is ) approximately inversely proportional to the square of ~ B D " ) . The average rate of reaction can be expressed in terms of ngD(f)as follows =
(9)
Again, if n B ~ ' f is ' assumed to be the only random variable, the variance of R is independent of n H D 1). In the hydrogenation studies under consideration here W / Q was l ) exheld nearly constant and different levels of n ~ ~ twere plored to obtain the kinetic data. Under these circumInd. Eng. Chem.,Fundarn., Vol. 13, No. 1, 1974
53
T a b l e 11. Fitting Models I and I1 to the Differential Data No. of data
Model
Temp,
I
O
points
K
365 375 365 375
I1
k
1 . 9 4 x 10-9 2 . 1 9 X 10-lo 1.36 X 2 . 1 5 X 10-l2
47 41 47 41
stances, E has an advantage of essentially constant variance as compared to the nonconstant variance of x i f ) . We shall outline the development of the regression model using rate model I1 as an illustration. The final form using rate model I is presented also. In an integral isothermal packed bed catalytic flow reactor, the steadystate fractional conversion of a reactant is given by the equation
The molar feed rate of butadiene can be expressed in terms of p ~ ~ ( 1 ) .
(11)
KBD
Statistic R2
0.96 0.19 0.0299 0.0232
0.982 0.929 0.961 0.929
Table 111. Integral Data. Results of Regression Analysis Best estimates Parameter ko KBD~ KBE A B
Residual error as fractional conversion of butadiene Statistic R2
Model I
Model I1
9 . 3 5 x 10-3 4 . 2 x 10-7 3 . 3 7 x 10-3 6958 - 4860 0 .0423
3.02 x 4 . 6 7 x 10-7 1 . 9 0 x 10-3 7460 - 4540 0.0421
0.9039
0.9049
T a b l e IV. Combined Data. Results of Regression Analysis Best estimates
We substitute Model I1 for R and note that for our integral data the following approximation h d d s PBD
= (1 -
X)PB,
(1)
(12)
For the present reactor system P H ~ ,PBE, and T are almost constant a t the inlet values. Therefore 22,400PW'=
'If'
+
KBJjpBE'l'
+
-
x)KBDPBD("13dr
k'(PH2('))l 5 p B D ' I ) ( l - x )
QPBd"
(13)
Parameter
Model I
Model I1
ko
1.11 x 10-8
KBD~ KBE
4.25 5.46
A
1923 - 4860 3 . 3 5 x 10-8
6 . 5 x 10-lO 1 . 4 7 x 10-7 -2.01 x 10-3 1599 - 4540 4.85 X
R (assumed) Residual error as g mole/ (sec) (g of catalyst) Statistic R 2
x 10-7 x 10-3
0.9631
0.9227
On integration we obtain
data. An additional advantage of treating the data by our method is that no line of demarcation between differential and integral data is required. However, if desired, for differential conditions, say -1% conversion in this case, the computations can be simplified since eq 8 is then a valid approximation.
where
Results of Regression Analysis
A similar treatment for model 1yields 22,400P W
-
QPBD(I'
(16) Equations 14 and 16 are implicit in the variable fractional conversion of butadiene. These equations can be solved by applying the Newton-Raphson root finding scheme. The obserued conversion serves as an excellent first guess. Since x t f icomputed from eq 14 or 16 is a predicted value, it will be denoted as k ( f J . The predicted value of the average rate of reaction, as defined by eq 7 , is then computed as 2
R=
Xf)pBD'I'&
22,400PW
(17)
Equation 17 coupled with eq 14 or 16 forms a regression model suitable for both integral and differential kinetic 54
Ind. Eng. Chem..
Fundam., Vol. 13, No. 1 , 1974
Since the integral data obtained by Kahney (1969) were not analyzed prior to this work, we analyzed them separately and then combined them with Gerheim's differential data (1966). An accurate estimate of the adsorption coefficient for the combined butenes was unnecessary because of its relatively small contribution to the denominator of each rate model. Therefore it was assumed to be constant, independent of temperature. The results of regression analysis for 48 integral data points are summarized in Table 111. The residual errors for both the models show clearly the disadvantage of using fractional conversion as a response variable for the combined data. The residual error for the integral data is about 0.04, which is four times as high as the actual conversion (-0.01) for the differential data. The regression models developed in the last section were used to analyze the combined data. This included all 145 data points reported in Table I. The results of the analysis are shown in Table IV. Four parameters, ko, KBno, K B E , and A , were estimated for both the models. Since the estimates of parameters A and B are highly correlated for the high butadiene partial pressure differential data and the separate estimates cannot be obtained directly, the estimate of B from analysis of the integral data was accepted as reliable.
Table V. Rate and Adsorption Coefficients from Various Types of Data
Data Model
Coefficient
I
I1
Temp,
Differential
O K
Integral
x 10-9 x 10-10
k k KBD KBD
3 65 375 365 375
1.94 2.19 0.96 0.19
k k KBD KBD
365 375 365 375
1.36 x 1O-I2 2 . 1 5 X 10-l2 0.0299 0.0232
4.92
Combined
x lo-"
5 . 7 2 X 10-11 6 . 5 8 X 10-11 0.258 0,181
lo-" lo-"
8 . 1 3 X 10-l2 9 . 1 4 x 10-12 0.037 0.027
8.18x 10-11 0.255 0.179 4.01 X 6.93 x 0.118 0.085
Table VI. Lack of Fit Tests"
Model I
__
Source of variation Residual Lack of fit (by difference) Experimental error a
S.S.
Model I1 M.S.
d.f.
1.579 X
141
1.12 X
1.5304 X
129
1.187 X
4.86
x
10-15
12
3.58
x
S.S.
d.f.
x
141
2.35
3.26 X
129
2.53 X
4.86 X
12
3.58 X
F 3.31
1O-lC
3.32
Critical value of F from statistical tables: F o 9~(120,12)= 3.4494; Fo.ls(= ,12)
It was possible to establish the temperature dependence within the models for the combined and integral data. With differential data separate analyses were performed by Kahney a t two temperatures. Table V compares the constants obtained from all the analyses. The values of rate and adsorption coefficients are reported a t each temperature so that comparisons can be made with the values obtained from differential data. For the preferred Model I, corresponding constants for combined and integral data are generally close ,in value, but those obtained from differential data differ greatly. Integral data were obtained a t low butadiene partial pressures where the surface of the catalyst is not saturated with butadiene; a t high partial pressure the surface is essentially saturated. Thus, the integral data are in the region which is critical for obtaining a correct model. Model Discrimination. Negative adsorption coefficients do not have any physical significance. Model I1 gives a statistically significant negative value for KBE. Therefore as a group the parameter estimates of Model I are more realistic than those of Model 11. Furthermore, the lack-of-fit test favors Model I. An estimate of error variance for all data was computed from the replicate runs. Its value is 3.58 x with 12 d.f. This means the estimated experimental error of a single observation is 1.89 x l o - * g mole/(sec) (g of catalyst). The analysis of variance tables for the lack of fit test for both the models are presented in Table VI. The null hypothesis that the lack of fit is zero cannot be rejected a t 99% level of confidence for Model I, while the hypothesis can be rejected for Model 11. On this basis Model I is accepted. Conclusions In this work we have presented a technique for analyzing corn bined integral and differential reaction kinetics data. The effectiveness of the technique has been demonstrated by applying it to the-rate modeling of butadiene hydrogenation. The variable R as defined in eq 7 has been shown superior to the fractional conversion when the combined data are being fitted. Extensive data on the hydrogenation of butadiene were available. A comprehensive rate model was developed for the reaction by combining the results of 97 differential and 48 integral rate data points summarized in Table I. The model over these conditions is given by the equation
=
10-13
3.3608.
.
M.S.
x
F
10-15
7.06*
This model should be considered more reliable than those fitted to different subsets of the data. The lack of fit test clearly favored this model over its rival. Because of the significant negative adsorption coefficient for butenes. Model I1 lacks physical meaning.
Nomenclature
A,B, C = parameters to account for the temperature dependencies of k , K B D , and KBE,respectively, "K F = molar feed rate of butadiene, g moles-sec k = rate constant, g moles/(sec) (g of catalyst) (mm)2.5 ko = preexponential factor of h k' = rate constant modified by the catalyst activity K B D = adsorption coefficient for butadiene, ( m m ) K B D= ~ preexponential factor of K B D , ( m m ) - l KBE = adsorption coefficient for butenes, (mm) K B E o = preexponential factor of K B E , (mm) - l M = a functional form of rate model ~ B D = number of moles of butadiene per mole of feed PBD = partial pressure of butadiene, mm p B E = partial pressure of butenes, mm P H * = partial pressure of hydrogen, mm P = total pressure in the reactor, mm Q = volumetric flow rate of the feed gases, cm3/sec at standard conditions R = rate of reaction, g moles/(sec) (g of catalyst) R2 = coefficient of multiple determination S = adsorption site on catalyst T = temperature, "K W = total weight of the catalyst in the reactor, g x = fractional conversion of butadiene y = 1-x. Greek Symbol N = defined in eq 15 Superscripts - = average value A = predictedvalue (f) = finalvalue (i) = initial or input value Ind. Eng. Chem., Fundam., Vol. 13, No. 1, 1974
55
Literature Cited Box, G. E. P., Hill, W. J., Technometrics 9,57 (1967). Froment, F. G., Mezaki, R., Chem. Eng. Sci. 25, 293 (1970). Gerheim, C. C., Ph.D. Thesis, Purdue University, 1966. Hill, W. J., Hunter, W. G . , Wichern, D. W., Technometrics I O , 145 (1968). Hougen,'O. A.. Watson, K. M., "Chemical Process Principles Part 1 1 1 , " p 906. Wiley, New York, N. Y.. 1943. Kahney, R. H., Ph.D. Thesis, Purdue University, 1969.
Kittrell, J. R.. Advan. Chem. Eng. 8, 97 (1970). Mezaki, R., Happel, J., Catal. Rev. 3, 241 (1969). Sane, P. P., Ph.D. Thesis, Purdue University, 1972. Sane, P. P., Eckert, R. E., Woods, J. M., Chem. Eng. Sci. 27, 1611 (1972).
Received for reczew February 5 , 1973 Accepted October 18, 1973
Axial Dispersion in the Turbulent Flow of Power-Law Fluids in Straight Tubes William B. Krantz" Department of Chemical Engineering, University of Colorado. Boulder, Colo. 80302
Darsh T. Wasan Department of Chemical Engineering, Illinois Institute of Technology, Chicago, Ill. 60616
The theory of Taylor for axial dispersion in turbulent flow in straight tubes is extended to power-law fluids. The effects of molecular and turbulent Schmidt number and soiute holdup in the wall region are incorporated in the analysis. The results agree well with axial dispersion data for the special case of Newtonian liquids. The predictions indicate that for Reynolds numbers larger than approximately 6000, pseudoplastic fluids will exhibit reduced axial dispersion, whereas dilatant fluids will exhibit enhanced axial dispersion relative to Newtonian fluids at the same flow conditions; for Reynolds numbers extending from the critical Reynolds number to 6000 the reverse trend is predicted. Presumably at larger Reynolds numbers the effect of the more blunt velocity profiles associated with decreasing flow behavior index predominates that of the thicker viscous sublayer, whereas at smaller turbulent Reynolds numbers the converse is true.
If a pulse of solute is injected into a fluid flowing in a straight tube it does not remain concentrated as a pulse but is axially dispersed relative to a reference frame moving at the average velocity. This dispersion is due to both axial diffusion and the velocity of the fluid particles relative to the average velocity. The former contribution is small compared to the latter for most fluids under a wide range of flow conditions. As the velocity profile approaches that corresponding to plug flow, the dispersion is reduced due to the decreased relative motion of the fluid particles. Radial diffusion also diminishes the axial dispersion since it reduces the radial concentration gradient. Interest in axial dispersion was generated by the need to measure the flow rates in long pipelines and to ascertain the degree of interdispersion when two liquids are transported successively in a pipeline. Flow rates in large water mains or in petroleum products pipelines can be determined by measuring the transit time of a tracer injected into the flow. Axial dispersion of the tracer yields an output concentration distributed in time; thus it is necessary to identify that point on the concentration-time curve which corresponds to the average velocity. Different liquids, for example regular and premium gasoline, are oftentimes transported in long pipelines successively without any separation. Intermixing due to axial dispersion produces a slug of liquid which must be purified or discarded in order to maintain product quality; it is of interest to predict the extent of the contaminated slug of material. Axial dispersion is also of importance in the design of tubular reactors since it decreases the driving force. The 56
Ind. Eng. Chem., Fundam., Vol. 13, No. 1, 1974
recent work of Gill and Sankarasubramanian (1971, 1972) on the dispersion of nonuniformly distributed time-variable sources in time-dependent laminar flow in straight tubes should permit many new applications of axial dispersion theory, a few of which are indicated by these authors. Axial dispersion in both laminar and turbulent flow of Newtonian fluids has been extensively studied. Sir Geoffrey Taylor (1953, 1954a) and later Aris (1956) developed the theory for axial dispersion in Newtonian fluids for laminar flow in straight tubes. Axial dispersion in turbulent flow was first analyzed by Taylor (1954b); however, his results are restricted to Reynolds numbers greater than 20,000 and to fluids having a Schmidt number of unity. Improved forms of the mean velocity profiles for both the turbulent core and the wall region which have been developed in recent years have allowed Taylor's analysis to be extended to the full range of turbulent Reynolds numbers. Tichacek, et al. (1957), extended Taylor's results to include large Schmidt numbers and used an improved form of the mean velocity profile based on averaging the velocity profiles of many investigators. Nonetheless their results may be inaccurate a t smaller turbulent Reynolds numbers since experimental velocity profiles are not available for the viscous sublayer in the region very close to the wall. Although this region is small in extent, it contributes greatly to the axial dispersion; thick viscous sublayers can hold up a considerable amount of solute, thus greatly increasing the axial dispersion. Flint and Eisenklam (1969)
