R E S P O N S E S U R F A C E M E T H O D S IN HETEROGENEOUS KINETIC MODELING J. R. K I T T R E L L Chevron Research Co., Richmond, Calif. 94802 J O H N ERJAVEC University sf Wisconsin, .Madison, W i s . 54537
In process development studies for reacting systems, a reaction model must usually be specified which adequately de!icribes the reaction kinetics; however, process improvement considerations often demand experiments different from those optimal for reaction rate modeling. A compromise can be obtained b y examining the empirical representations of the important responses, obtained via response surface methodology. These empiirical representations of the response surfaces can directly assist in the specification of a reaction rate model,, indicate regions of further experimentation desirable in the modeling program, and provide information necessary for process improvement and preliminary process design and economic considerations. The discussion i s illustrated by data on the vapor-phase isomerization of n-pentane in the presence of hydrogen and over a supported metal catalyst.
HE determination of a model which adequately represents T a reacting system has been the object of considerable effort in heterogeneous kinetics. Numerous techniques have been proposed, depending upon the type of model used (Kittrell et al., 1966a). Every modeling technique, however, must examine the nature of the response surface under consideration. For example, the approach of Yang and Hougen (1950) requires an examination of the dependence of the initial reaction rate upon total pressure; a maximum in this cross section of the reaction rate surface may suggest that only one or a few dual-site surface reaction-controlled mlodels need be considered further. This paper indicates the advantages of using response surface methods rather than other existing modeling procedures for an inchoate study of the plausible kinetic models. T h e primary advantage of response surface methodology is that it allows an adequate fit and interpretation of the multidimensional response surface, instead of, for example, only the unidimensional pressure-dependent data mentioned above. This multidimensional surface representation specifies a functional form which an adequate reaction rate model must possess, regions of experimentation necessary for further modeling and process improvement studies, and an efficient empirical model for preliminary process development activities.
Existing Modeling Procedures
One class of models used for fitting heterogeneous kinetic data, advocated by U’eller (1956), is the power function model: =
kiPAnPBrn
(1)
The methods generally used for analyzing data with such a model are largely those used for years by chemical kineticists (Kittrell et al., 1966b). These methods are developed primarily for a single reacting component, for which Equation 1 may be written as
r = kpAn
(2)
If combined with the equation of continuity for the reactor and integrated, the adequacy of assumed reaction orders may be tested by the method of integration. Alternatively, reaction
rates may be obtained and the method of differentiation used to test the adequacy of assumed reaction orders. I n either case, one often seeks to rearrange the model into a linear form and plots the data to ascertain their ability to form a linear graph. Other methods have employed fractional lifetimes, dimensionless curves, and statistical transformations (Kittrell et al., 1966b). These and other numerical techniques, such as linear least squares, are similar in that they effect an algebraic transformation of variables so that the shape of the appropriate response surface can be interpreted in an unambiguous fashion. For these linearizable one-variable cases, response surface methodology can add very little to the experimenter’s ability to analyze the data. For the cases in which several reacting components are present (Equation l ) , the procedure is not so straightforward. One method of data analysis in such situations requires the measurement of reaction rate data at such large concentrations of all except one reactant that their concentrations effectively remain constant during the entire course of the reaction. Then, the order of the low concentration reactant can be determined by any of the single-component methods. This isolation method allows a determination of the component reaction orders, but a very limited region of the experimental space has been covered in determining these orders. Since the model has not been tested for experimental conditions in which all concentrations are simultaneously varying, it should be used with caution here. These estimated constants, on the other hand, should provide good initial estimates for a nonlinear least squares program which can analyze data taken when the concentrations of all the components are varying. Initial reaction rates may also be used for fitting data with Equation 1. Here, the initial concentrations of the individual reactants are varied while all other reactant concentrations are held constant. T h e data can then be analyzed by the singlecomponent methods. An advantage of this technique over the isolation method is that the concentrations of the reactants can be nearly equal instead of some concentrations being in large excess. The method allows reaction rates to be obtained over the entire range of composition with respect to known VOL. 7
NO. 3
JULY 1 9 6 8
321
EQUATION 3
?\
X-
W
m
5
a
a
> I2 W
a
z
P
W
W e
P
t /
TOTAL PRESSURE,
; / x2
?f
\
INDEPENDENT VARIABLE x2
Figure 1 . Hypothetical initial rate dependence upon total pressure
major reaction participants. However, these rates may not be equal to those obtained from experiments in which extensive conversions are allowed, because of the presence of trace byproducts during the reaction which affect the reaction rate. (Such considerations have prompted the definitions of “order with respect to concentration” and “order with respect to time.”) A more commonly published method of analysis in heterogeneous kinetics uses Hougen-Watson or LangmuirHinshelwood models. In such cases we typically attempt to discriminate among large classes of possible models (such as adsorption or surface reaction-controlling models) or among several surface reaction-controlling models, such as the single site
r=-
kKs(xz Kixi
(1
+
(1
+ Kixi +
+
K2~2
(3) K3~3)
or the dual site
r =
kK2(x2 - x 3 / K )
+K3d2
(4)
K2~2
Here, reaction data are often taken as a function of the component partial pressures, XI, X Z , and x 3 . From these data, the initial reaction rate is extracted and plotted as a function of total pressure, such as shown in Figure 1 where the initial rate data suggest that the dual-site model is preferred. If several models can describe the initial rate data of Figure 1, as would be the case if data were taken only up to pressure rl, these models are linearized by an algebraic transformation and fitted to all of the available data by linear least squares. Models possessing estimates of the adsorption constants which are definitely negative are eliminated on physical grounds. This entire procedure has been discussed in detail (Kittrell et al., 1965). One major disadvantage of this method, the isolation method, and the method of initial reaction rates is the necessity of using a onefactor-at-a-time experimental design to allow the several twodimensional graphs to be drawn, Response surface methods, by revealing the major features of the reaction rate response surface, can provide this same information without such a restriction. More sophisticated methods of selecting experiments for 322
Contours of constant reaction rate
discriminating among several rival models have been proposed (Box and Hill, 1967; Kittrell et al., 1966a). The number of possible Hougen-Watson models taken into consideration in selecting these experimental designs should, practically speaking, be less than about 10 to apply such procedures effectively. However, more than 80 such models have been postulated for methane oxidation (Hunter and Mezaki, 1964). Hence, experiments must be conducted to remove the obviously inadequate models from consideration early in the modeling program. By allowing an interpretation of the reaction rate response when all variables are changed simultaneously, inadequate models can be eliminated through response surface methods. Response Surface Methodology
- x3/K)
+
Figure 2.
l & E C PROCESS D E S I G N A N D DEVELOPMENT
Two primary features of response surface methodology are the experimental design and the method of analysis (Box, 1954; Davies, 1960). I t is advantageous to use one of the common first-order or second-order designs (such as the factorial or composite designs) with response surface methods. I n kinetic applications of response surface methodology, an empirical equation is desired which will approximately describe the kinetic rate surface. One conceivable surface, for two independent variables, could have the contours of constant rate as shown in Figure 2. If the nonlinear reaction rate model can be written
r
=
f(x1, X P ;
K)
(5)
the surface of Figure 2 can often be represented approximately, in the region where data are taken, by the Taylor expansion
x = xs
x = xs
x = XQ
x = XI
or as 7
+ bixi +
= bo
+
b 2 ~ 2
+
612x1~2
b1ixi2
+
622x2'
(7)
This is the general equation for an ellipse in two independent variables; in order for Equation 6 to approximate Equation 5, the contours of Figure 2 must be approximately elliptical in the region of experimentation. If a new coordinate system may be found which is centered about point S and is rotated to the directions denoted by axes X I and X2, then Equation 7 simplifies to 7
- rS
=
BiiXx2
+ B22Xz2
(8)
where rs is the rate at the stationary point, S. This transformation of axes, a canonical transformation, has eliminated the crossproduct terms of Equation 7 . From Equation 8, then, the nature of the surface can be ascertained by examining the signs and magnitudes of coefficients B11 and B22. The 10-term equation for three independent variables, corresponding to Equation 7 , reduces to
For example, if B11, B22, and B33 were all significantly positive, any movement along the independent variables X i should result in an increase in rate 7. Thus, the surface has a minimum a t point S. Similarly, if all the B's are significantly negative, a maximum exists at S. A saddle point exists if the B's are of unequal signs, while iS some are zero a ridge system exists. A more complete description of such matters has been presented (Box, 1954; Davies, 1'160). T h e mathematical (details of carrying out a canonical transformation have also been presented (Box, 1954; Davies, 1960). For the present purposes, one must first find the stationary point, S-e.g., by differentiation of the fitted Equation 7. Then the rotation of the axes is achieved by finding the latent roots and vectors of a matrix of the estimated coefficients of the quadratic and interaction terms of Equation 7. The latent roots are, in fact, the coefficients B i t ; and the latent vectors relate the rotated to the unrotated axes. This relation is usually summarized in the forin: XI
-
XlS
This means that any axis axes XI, x2, and x 3 by
xi
=
q i l h
-
XlS)
Xiis related to the untransformed
+ q.12bz -
22s)
+
qis(X3
-
x3s)
(10)
and conversely that xi
- xis
:= q1iX1
+
q2&2
+
q3iX3
(11)
The requirement that these two sets of axes be linearly related is obvious from the rotation and shift of the axes shown in Figure 2. Response Surface Methods and Kinetic Modeling
A detailed example of the use of the response surface knowledge to gain theoretical insight into a kinetic equation has been published (Box and Youle, 1955). The use of such methods to investigate the influence of several process variables and catalysts on reactions has been discussed (Franklin at al., 1956, 1958). Power Function Models. One obvious application of response surface techniques to models such as Equation 1 is to indicate the gross adequacy of the model to represent the rate
data. I n particular, the model of Equation 1 cannot describe a surface that exhibits a finite stationary point along axes pa, p B , etc. With response surface methodology, the exhibition of a stationary point in a multidimensional reaction rate surface is immediately apparent. Furthermore, if the data taken u p to the time of the analysis do not exhibit a stationary point but the fitted surface indicates that a maximum may exist slightly outside the experimental region, data could be taken in this latter area to determine whether the stationary point is a real feature of the reaction mechanism or only an erroneous implication of the empirical fit of the data. I n this fashion, not only could the range of applicability of Equation 1 be roughly outlined, but also the existence of regions where experiments should be conducted to indicate the true nature of the surface could be suggested. If a stationary point exists, the model of Equation 1 must be changed. Response surface techniques can also be used to determine the order of a reaction. This has been discussed in detail by Pinchbeck (1957). Hougen-Watson Models. A very fruitful area of application of response surface methods is the preliminary screening of rival Hougen-Watson models. A common method of conducting. this initial screening is by inspection of the initial rate-total pressure profile; it allows a visual inspection of a cross section of the entire reaction rate surface which can be directly related to various model types. There is no reason to be restricted to a n inspection of this cross section of the surface, however, because response surface methods allow an examination of the entire surface utilizing, say, a central composite design rather than this one-factor design (Kittrell et al., 1965). Instead of examining each model for the ability of its initial rate us. pressure to exhibit a maximum, for example, the model could be tested for its ability to exhibit a stationary point along the axis of every independent variable. For example, Equation 3 does not exhibit a stationary point along any axis x i except at equilibrium, whereas Equation 4 should exhibit a saddle point. Additionally, the improved coverage of the reaction rate surface by the composite design will generally allow a more effective data analysis than that obtained by any of the one-factor-at-a-time experimental designs (Kittrell et al., 1965). One primary advantage of the application of response surface methodology to kinetic modeling, then, is that it allows an interpretation of the entire reaction rate surface in the modeling program-not just one section of this surface obtained by a one-variable-at-a-time design. Not only can this suggest a preferable mechanistic model, but it can also indicate regions of experimentation valuable in the future modeling program. A very important second advantage is also apparent. Experimental programs for mechanistic modeling can be lengthy in terms of the total time elapsed before a satisfactory conclusion is reached. I n many process development studies, however, the preliminary process design and economic evaluations cannot wait until the modeling program is completed. I n fact, it would be unwise to invest heavily in a mechanistic modeling program before the preliminary process design and economic evaluations indicate that the over-all process has an acceptable probability of being economically successful. Being statistically sound, the use of response surface methodology early in the modeling program will provide one of the best empirical models that can be obtained. I n addition to representing the existing data well, it will indicate regions of high yields that should be experimentally examined in further process development studies before the mechanistic modeling program is begun in earnest. Hence, not only does response surface methodology VOL. 7
NO. 3
JULY 1 3 6 8
323
br = -0.49 Table I. 21
R~~ No. 1
x1
=
- 300 100
-1 $1
2
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
Experimental Design
-1 +1 -1
+I
-1 $1
-1.68 $1.68 0 0 0 0 0
32
~2
-
50 -1 -1
1 1 -1 -1
3s =
- 84
50 -1 -1 -1 -1
1 1 0 0
0
-1.68 $1.68 0 0
0
0
0
0 -1.93 1.17 -0.49 -0.50 -1.55
%a
1 1 1 1 0
0 0 0
0
=
150
-1.31 1.75 -1.29 3.56 0.18 3.68
0 0 -1.68 $1.68 0 0 0 0 -0.38 -1.57 -1.07 -1.53 -1.61 -1 3 2
Example
Data have been published (Carr, 1960) on the isomerization of n-pentane to 2-methylbutane, taken according to an augmented central composite design. These were then analyzed by linear least squares to find an appropriate Hougen-Watson model. We will illustrate the use of response surface methods to achieve the same end while simultaneously obtaining an empirical model useful for preliminary process design studies and determining areas of fruitful future experimentation. T h e reaction rate data were taken with a flow differential reactor. The only reaction considered was the isomerization of n-pentane to 2-methylbutane, since the equilibrium conversion of neopentane is small at the reaction temperature of 750' F. Hydrogen was present as an inert, resulting in three independent variables. The catalyst and the reactor have been described (Carr, 1960); the coded data are reproduced in Table I. Now consider an empirical representation of the rate surface to ascertain what general characteristics a good model should possess. One empirical model fitted to the rate data for this purpose is: =
bo
+ 6121 + 6222 + b323 4- biiZi2 + b2222' +
63333'
+
-3% --
Subscripts 1, 2, and 3 refer to hydrogen, n-pentane, and 2methylbutane, respectively. The parameter estimates and their standard errors, obtained by a linear least squares fit of Equation 12 to the data, are shown in Table 11. The analysis of variance of Table I11 suggests that this equation, including the second-order terms, adequately represents the data. A general interpretation of the shape of the rate surface may be obtained from this equation. To do this, we can differentiate Equation 12 which, using the coefficients of Table 11, yields l&EC PROCESS DESIGN A N D DEVELOPMENT
-0.154
bj.12
br = 1.26
- 0.241 21 - 0.416 Zz - 0.158 Zs
bj.2
bqr
= -0.416
2
-b- -
-1.50
bj.8
+ 0.505 - 0.158 + 0.514 3s 21
--
bj.2
-
$2
0.514
Hence, a stationary point exists at 21 = 1.32, 22 = 1.43, and 2.08, found by a simultaneous solution of Equations 13, 15, and 17 set to zero. An examination of Equations 13 through 18 indicates that the surface exhibits a saddle point, with a maximum along the 2 2 axis and a minimum along the 23 axis. I t is believed that Equations 13 and 14 merely reflect a zero rate at equilibrium, since if 6 1 1 = 0 the equations simplify nearly to an equilibrium relationship. The importance of the precise estimation of the coefficients of the interaction terms, b i j , can be seen from the role they play in estimating the location of the saddle point. These interaction terms cannot be estimated from one-factor-at-a-time experimental designs. One cross section of the saddle point predicted by Equation 12 is plotted in Figure 3. Given that the rate data suggest the existence of a saddle point, let us consider the models of Table IV, derived using the formalism of Yang and Hougen (1950). These are not all of the possible Hougen-Watson models which could describe this j.3 =
2.c
1,75
1 2
w" a v) 3 v) W
1.5C W 2
3z W
4 n
8
324
+ 0.505 28
$2
bz 1 Reaction Rate, Grams/ Gram Cat.-Hr. 3.541 2.397 6.694 4.722 0.593 0.268 2.797 2.451 3.196 2.021 0.896 5.084 5.686 1.193 2.648 3.303 3.054 3.302 1.271 11.648 2.002 9.604 7.754 11.590
provide a n effective way of approaching a mechanistic model; it also has numerous engineering advantages important from the process design and economic standpoints.
Y
- 0.154 21 - 0.241
1.25
I.o
I .75
-
2 ,o
2,25
CODED i PENTANE PRESSURE, 13
Figure 3.
Predicted saddle point of reaction rate Coded hydrogen pressure, XI = 1.32
Table II.
Parameter Estimates for Fitted Models
Estimated Value Equation 79 -2.51 zk 10.8 0.0676 =t 0.024 0.00519 f 0.037 0.185 f 0.052
Equation 12 3.06 0.20 0.12 -0.490 1.26 Z!C 0.12 -1.50 f 0 . 1 2 -0.077 f 0.12 -0.208 f 0.07 0.257 0.12 -0.241 f 0.12 0.505 0.14 -0.158 =t 0.14
Parameter
+ + + +
baa biz bia b23
Table 111.
Source Equation 12 Mean squares Degrees of freedom
Extra Due to First-Order Terms
76.4 3
Analysis of Variance for Rate Equations
Extra Due to Second-Order Terms
Residual
1.72 6
0,225 14
c
Mean square ratios F statistics 95 % 97.5% Equation 19 Mean squares Degrees of freedom
-
d
L
I
d
Y
Adsorption-controlled with dissociation
d
3.1 3.9
8.7 14.3
Model k(Xz - x s / K )
[1 + Kixi ($ +
[
k(xz
+ Klxl +
K3) X I ]
- x3/K) -
q< + x3
Rate Surface for Model No stationary point except at equilibrium; first derivative with x i zero only at equilibrium; first derivative with xp always positive and x3 always negative except at equilibrium
- x3/K)
kKz(xn
-x~/K)
No stationary point except at equilibrium; first derivative with x i zero only at equilibrium; first derivative with xp always positive and x t always negative except at equilibrium Can have saddle point; first derivative with x i zero at equilibrium; first derivative with x z may be zero with second derivative negative; first derivative with x3 may be zero with second derivative positive Same as 4
kKz(x2
- xa/K)
Same as 4
4. Dual site, surface reaction-controlled
=
('
+ K1xl + K2X2 + K3x3)2
=
(1
+ d K X + K Z X +Z Ksxs)'
kKs(xz
+ Kixi + K Z X Z+ kKz(Xz
7
d
No stationary point except at equilibrium: first derivative with X I zero only at equilibrium; first derivative with x 3 may be zero with second derivative positive K8xa]
(1
7. Desorption-controlled
Y
L
5.9
=
K3~3)
-x~/K)
+ Kixi + d K X + Ksxal' kK(x2 - x s / K ) = [l + Kixi + ( K z + KKa)xzl
r =
0.090 3
13.0
3. Single site, surface reactioncontrolled
Surface reaction-controlled dissociated n-pentane
0,528 17
Possible Reaction Rate Models for Pentane Isomerization
1
5. Dual site, surface ,reactioncontrolled, hydrogen dissociated
d
8.7 14.3
Y
r =
L
3.1 3.9 0,462 20
1. Adsorption-controlled single site r =
10.88 3 10.6
6,019 3
IV.
115.56 117
8.7
L
k
8.8 14.3 99.86 20
Table
0.0954 3 2.73
868.41 3 L
Pure Error
0.260 11
2.9 3.5
Mean square ratios F statistics 95 % 97.5%
6.
Lack of Fit
7.66
Mean sauare ratios F statisti'cs 95 % 97,570 Equation 20 Mean squares Degrees of freedom
2.
Equation 20 2.58 =!= 0.74 0.00577 f 0.0016 0.0009 =t0.0025 0.0153 i 0.0035
[l
system, and not all of the models of Table IV are equally plausible. For example, the dissociation of pentane is unlikely and the adsorption or desorption controlling models would not be expected to be appropriate at the relatively low pressures of this study. The descriptions of the rate surfaces of Table IV suggest that only the dual-site surface reaction models predict the proper saddle point, the simplest of which is model 4.
Same as 1
Hence, using the knowledge of the rate response surface alone, the data would tend to favor the dual-site model of Equation 4, since it is the simplest model predicting a saddle point. However, because of the inadequate coverage of the stationary point by the experimental data (note the position of the stationary point relative to Table I), it cannot be said with any great degree of certainty that the model of Equation 4 is VOL. 7
NO. 3
JULY 1968
325
Table V.
Equation 3
Equation Form Nonlinear
4
Nonlinear
19 20
Results of linear and Nonlinear least Squares Analysis for Single and Dual Site Models
k
K1
36.3 2c 7.7 133.0 =k 20
0.0619 f 0.14 0.0022 =k 0.00037 -0.0269 0.00222
193 426
Linear Linear
Kz 0.0327 f 0,076 0.0013 =k 0.00037 -0,00207 0.00035
preferred. Apparently, additional data in the region of the predicted stationary point are required to pinpoint the model precisely. Hence, we are provided with information as to the experimental design needed in the next step of the modeling procedure. The use of more quantitative experimental design procedures (Box and Hill, 1967) will also suggest that these experiments be conducted first. Also, we have an indication of the better mechanistic model of those considered and a good empirical representation of the surface for process development through Equation 12. Before considering further the information to be gained from the empirical surface representation, let us compare the tentative conclusion of the preference for model 4 of Table I V with the conclusions that could be obtained from a detailed analysis of the most plausible models of Table IV, models 3 and 4. The inclusion of models 5 and 6 in the analysis will not change the conclusions below. The models can be linearized by an algebraic transformation to yield :
K3
0.146 f 0.32 0.0049 f
Ratio F to Pure Residual Error Statistic M e a n Square Mean Square (95%) 0.162 1.69 8.7
0.174
1.81
8.7
99.9 0,462
9.18 5.13
8.7 8.7
0,00089
-0.0737 0,0059
I t is reasonable to expect that data in the region of the predicted stationary point would allow the rejection of one of these models. At this stage, a great deal has been learned about the modeling requirements for the pentane isomerization through response surface methods. First, the dual-site model 4 can be used as an adequate model, although extrapolations beyond the region of experimentation are particularly dangerous. The same conclusion has been drawn by a more detailed modeling analysis. I n addition, a consideration of the response surface has suggested regions of experimentation (near the stationary point) which should be examined if it is desired to pinpoint the model more precisely. Additional information for process improvement can be obtained by writing Equation 12 in its canonical form: r
-
2.11 = 0.155 Xi2 - 0.292 XzZ
+ 0.419 Xa2
(21)
The standard errors of the coefficients of X?,Xz2, and Xa2 are, respectively, 0.094, 0.108, and 0.095. The canonical variables are related to the untransformed variables by : R-I x 1
xz x 3
- 7.32 0.599 0.639 0,483
R-y
-
- 2.08 -0.489
1.48
R-a
-0.635 0.746
-0.187 0,852
-0.200
For example, -xz+-
.\/kKz =
bo
K 3 x3
4%
(20)
+ 61x1 4- bzxz + b 3 ~ 3
These models were fitted to the data by unweighted linear least squares, resulting in the parameter estimates of Tables I1 and V. The unweighted analysis is justified only in that no better weighting estimate is available and the variation in reaction rates in Table I is relatively small (Mezaki and Kittrell, 1967). The analyses of variance of Table I11 indicate that both models adequately fit the data, Equation 20 somewhat better than Equation 19. The negative estimate of the bo term of Equation 19 in Table I1 (and, hence, the negative estimate of Ki in Table V) is of no great importance, since the confidence interval is relatively large. Equations 3 and 4 were also fitted to the data by unweighted nonlinear least squares, resulting in the parameter estimates shown in Table V. The residual mean squares, when compared to the pure (experimental) error mean square, again indicate that it would be dangerous to reject either model. No unusually severe residual trends exist for either model when fitted by nonlinear least squares. All of the parameter estimates are now positive. The preferred parameter estimates for these two models are collected under the nonlinear equation forms of Table V. The conclusions of this detailed analysis are thus identical to those obtained by response surface methods. 326
I&EC PROCESS DESIGN A N D DEVELOPMENT
Xi
= 0.599 (31 - 1.32)
-
0.635 (32 - 1.48) 0.489
-
(33
-
2.08)
(22)
From Equation 21 and the associated standard errors, it can be seen that the variable XIaffects the rate slightly, if at all. A saddle point in the variables XZand X3 is also indicated by this equation. This indicates that, for process development work, the variable X3 should be explored to increase the performance of the process. If some response other than reaction rate is considered to be more indicative of process performance, such as cost, yield, or selectivity, the canonical analysis corresponding to Equation 21 could be performed on this response to indicate regions of the independent variables to be explored for increasing process performance. Since these responses of cost, yield, and selectivity are generally calculable from the data obtained from the design of Table I , this process development insight is gained in addition to the mechanistic modeling considerations just discussed. The utility of response surface methods for such process development studies is widely recognized (Box, 1954; Davies, 1960). Hence, when data are being taken for the modeling of a chemical reaction in process development studies, it is extremely valuable to begin the experimentation with an examination of the response surfaces of interest, whether a rate surface as considered here or a yield or selectivity surface. If the response surface analysis were done by hand, so much time would
= reaction rate, grams/gram catalyst-hour = reaction rate at stationary point, grams/gram catalysthour = canonical axis, i, i = 1, 2, 3 vector of all of variables, Xi = partial pressure of hydrogen, atm. = partial pressure of n-pentane, atm. = partial pressure of 2-methylbutane, atm. = magnitude of variable x i at stationary point, i = 1, 2, 3 = coded values of variables x i defined in Table I, i = 1, 2, 3 = vector of components, x i a , i = 1, 2, 3 = total pressure, atm. = pressure level of Figure 1, atm.
be required that the attention of the development engineer could be taken from the other equally important engineering matters. However, now "canned" response surface programs exist which can be used with little time investment by the engineer. O n a digital computer, for example, the complete response surface analysis and a partial printing of contours of the response surface (such as that shown in Figure 3) were completed in about 1.5 minutes. Hence, response surface methods would not be a diversion, but rather an important tool to interact with and to guide the chemical engineering judgment in a modeling program within a process development study.
r rs
Nomenclature
literature Cited
B,,
Box, G. E. P., Biometrics 10,No. 1,16 (1954). Box, G. E. P., Hill, l V . J., Technometrim 9 No. 1 , 57 (1967). Box, G. E. P., Youle, P. V., Biometrics 11, No. 3, 287 (1955). Carr, N. L., Ind. Eng. Chem. 52, 391 (1960). Davies, 0. L., "Design and Analysis of Industrial Experiments," 2nd ed., Hafner Publishing Co., New York, 1960. Franklin, N. L., Pinchbeck, P. H., Popper, F., Trans. Inst. Chem. Em. 34. 280 (1956). Franglin, 'N. L.', Pinchbeck, P. H., Popper, F., Trans. Inst. Chem. Eng. 36, 259 (1958). Hunter, l V . G., Mezaki, R., A.I.Ch.E. J . 10,315 (1964). Kittrell, J. R., Hunter, W. G., Watson, C. C., A.I.CI2.E. J . 11, 1051 (1965). Kittrell, J. R., Mezaki, R., lb'atson, C. C., Brit. Chem. Eng. 11, No. 1. 15 11966a). Kittrell,' J. R., Mdzaki, R., \\'atson, C. C., Ind. Eng. Chem. 58, S o . 5, 50 (1966b). Mezaki, R., Kittrell, J. R., Ind. Rng. Chem. 59, No. 5, 63 (1967). Pinchbeck, P. H., Chern. Ene. Sci. 6,105 (1957). T.\'eller, S.,'A.I.Ch.E. J . 2, 53 (1956). Yang, K. H., Hougen, 0. A, Chem. Eng. Progr. 46, No. 3, 146 (1950).
canonical parameters defining nature of the response surface in transformed coordinates, i = 1, 2, 3 6, = parameter values defining linear effects of the response surface, z = 0, 1, 2, 3 6,, = parameter va1ue.s defining the quadratic and interaction effects of the response surface, i = 1, 2, 3;. j = 1, 2, 3 f = abbreviation of the complete modelf(x; K) in Equation 6 thermodynamic equilibrium constant for reaction of npentane to 2-1nethylbutane, equal to 1.632 equilibrium adsorption constant for hydrogen, atm.-' equilibrium adsorption constant for n-pentane, atm.-' equilibrium adsorption constant for 2-methylbutane, atm.-l vector of all parameters in a model to be estimated from kinetic data forward rate constant, grams/gram catalyst-hour empirical constant in Equations 1 and 2 reaction order with respect to any component A reaction order with respect to any component B partial pressure of any component A partial pressure of any component B components of latent vectors, relating transformed coordinate system and untransformed system =
xi x = XI x2 x3 Xia
Ri x, 7 r 7r1
RECEIVED for review June 17, 1966 RESUBMITTED December 28, 1968 ACCEPTED March 16, 1968
GAS-SOLID HEAT TRANSER IN FLUIDIZED BEDS R. S . M A N N A N D L . C. L . F E N G Department of Chemical Engineering, University of Ottawa, Ottawa, Canada
The unsteady-state heat transfer between gas and solid particles was studied in fluidized beds, 2 and 4 inches in diameter. A system of transient heating and cooling of glass beads, silica gel, and alumina between 130" and 288" F. was used. The effect of several variables, particle size (0.004to 0.241l inch), bed settled heights (0.8to 16 inches), particle densities (8.5 to 206 pounds per cu. foot), thermal conductivities [0.013to 1.8 B.t.u./(hr.)(sq. ft.)(' F./ft.)], and air velocities (0.543to 4.347feet per second) on the spaceaveraged heat transfer coefficient, U, was investigated. Using a new approach in interpreting the driving force, a correlation for two ranges of Reynolds numbers, 10 to 60 and 60 to 2200, was developed.
NE
of the inherent properties of fluidized beds is the high
0 rate of heat transfer between the beds and the fluidizing medium. The transfer of heat between particles, fluids, and the surfaces in contact with them is as complex in its facets, and even more so in its mechanisms, as the problems associated with the many phases of fluid flow in such two-phase systems. T h e phenomenon of heat transport in fluidized beds has been the subject of numerous studies (Zenz and Othmer, 1960). Basic equations for fluidization (Leva, 1959), heat transfer from bed to wall and vice versa (Heerden et al., 1951 ; Leva and FVeintraub 1949; Toomey and Johnstone, 1953), and heat transfer in pneumatic fluidized systems (Koble et al., 1951;
Richardson and Ayers, 1959) were comparatively well investigated. The single-par ticle technique (Johnstone et al., 1941), theoretical considerations (Zenz and Othmer, 1960), and the early investigations (Heertjes and McKibbins, 1956; Kettenring et al., 1950; Shakhova and Rychkov, 1957; Walton et al., 1952) did not provide enough information to understand the complex phenomenon of heat transfer in fluid-solid systems. More recent works (Ferron, 1961 ; Frantz, 1961; Fritz, 1956) have brought a new approach to the problem of gas-solid heat transfer in fluidized beds. However, the work of Frantz was limitvd to liquid fluidizing medium and that of Ferron and Fritz to one solid sample only. VOL. 7
NO. 3 J U L Y 1 9 6 8
327
