I n d . Eng. Chem. Res. 1987,26, 2092-2098
2092 I.(
Laddha, G. S.; Degaleesan, T. E. Transport Phenomena in Liquid Extraction; McGraw-Hill: New York, 1978. Lahiere, R. J. Ph.D. Dissertation, The University of Texas at Austin,
= viscosity, CP (0.01 g/(cms))
constant, 3.1416 p = density, g/cm3 CJ = interfacial tension, dyn/cm 4 = operational holdup, dimensionless r =
1986.
Subscripts
c = continuous phase
d = dispersed phase f = formation r = rise (or fall) s = superficial cross-section basis S y s t e m Identification A = acetone
E = ethanol I = 2-propanol LC = carbon dioxide, near-critical or liquid SC = carbon dioxide, supercritical T = toluene W = water Registry
No. COz, 124-38-9; H,CCH,OH,
64-17-5; H 3 C C H -
(OH)CH,, 67-63-0.
Lewis, W. K., Jr. Ind. Chem. Eng. 1936, 28, 399. Major, C. J.; Hertzog, R. R. Chem. Eng. Prog. 1955, 51, 17. McHugh, M. A.; Krukonis, V. J. Supercritical Fluid Extraction, Principles and Practice; Butterworth: Stoneham, MA, 1986. Moses, J. M.; de Filippi, R. P. "Critical-Fluid Extraction of Organics from Water" Contract DE-AC01-79CS40258, June 1984; US. Department of Commerce, Springfield, VA. Paulaitis, M. E.; Krukonis, V. J.; Kurnik, R. T.; Reid, R. C. Reu. Chem. Eng. 1983, 1, 179. Peter, S. Ber. Bunsenges. Phys. Chem. 1984,88, 875. Pilhofer, T. Chem. Eng. Commun. 1981, 1I , 241. Pilhofer, T.; Goedl, R. Chem.-Znq.-Technol. 1977, 49, 431. Popovich, A. T.; Jervis, R. E.; Trass, 0. Chem. Eng. Sci. 1964, 19, 357. Rathkamp, P. J. M.S. Thesis, The University of Texas at Austin, 1986.
Rocha, J. A. Ph.D. Dissertation, The University of Texas at Austin, 1984.
Rocha, J. A.; Humphrey, J. L.; Fair, J. R. Ind. Eng. Chem. Process Des. Deu. 1986, 25, 862. Seibert, A. F. Ph.D. Dissertation, The University of Texas at Austin, 1986.
Literature Cited Brunner, G., Presented at San Francisco AIChE Meeting, Nov 1984. Brunner, G.; Kreim, M. K. Ger. Chem. Eng. 1986, 9, 246. Eisenbach, W. 0. Proc. Am. Chem. SOC.1985, 30(3), 206. Gilbert, M. I,. M.S. Thesis, The University of Delaware, Newark, 1985.
Handlos, A. E.; Baron, T. AICHE J . 1957, 3, 127. Higbie, R. Trans. AZChE t935, 31, 365. Hough, E. W.; Heuer, G. J.; Walker, J. W. Pet. Trans. AZME 1959, 216, 469. Johnston, K. P. In Encyclopedia of Chem. Tech., 3rd Ed.; Wiley: New York, 1984; Vol. 3, p 872. Kander, R. G. Ph.D. Dissertation, The University of Delaware, Newark, 1986.
Skelland, A, H. P.; Conger, W. L. Znd. Eng. Chem. Process Des. Deu. 1973, 12, 448. Thornton, J. D. Chem. Eng. Sci. 1956, 5, 201. Tiegs, C. Ph.D. Dissertation, The University of Erlangen, W. Germany, 1986. Treybal, R. E. Liquid Extraction, 2nd ed.; McGraw-Hill: New York, 1963.
Treybal, R. E. Mass Transfer Operations, 3rd ed.; McGraw-Hill: New York, 1980. Wilke, C. R.; Chang, P. AZChE J . 1955, 1 , 264. Received for review February 9, 1987 Revised manuscript received June 30, 1987 Accepted July 23, 1987
Use of Bifurcation Map for Kinetic Parameter Estimation. 1. Ethane Oxidation Michael P. Harold* and Dan Luss University of Houston, Department of Chemical Engineering, Houston, Texas 77004
Kinetic and transport parameters may be estimated from a set of experimental ignition and extinction points, i.e., a bifurcation map. Initial estimates of some of the parameters can be obtained from a preliminary analysis of temperature rise on the extinguished branch. The technique was used to fit the bifurcation points for ethane oxidation on a single Pt/A1,0, pellet. T h e dependence of the gas density on ambient temperature was found to have an important impact on the multiplicity features. T h e fit of data removed from the ignition or extinction points was not as good as that of the singular points. A method is presented for overcoming this limitation. It is well-known that steady-state multiplicity may occur in chemically reacting systems. The qualitative features of the multiplicity patterns may be used to infer the functional forms of the rate expression (Harold et al., 1987a,b). Moreover, the ignition and extinction points are sensitive functions of the kinetic parameters and can be used for kinetic model identification and fit. Hiam et al. (1968),Schwartz et al. (1971),and &der and Weller (1974) determined the activation energy of several reactions from measurements of the ignition points. Jiracek et al. (1971)
* Present address: Department of Chemical Engineering, University of Massachusetts, Amherst, MA 01003. 0888-5885/87/2626-20!32$01.50/0
used bifurcation diagrams, exhibiting multiplicity, to check the validity of a rate expression for hydrogen oxidation. Beskov et al. (1976) used ignition and extinction conditions to fit kinetic parameters for ammonia oxidation. Recently, Herskowitz and Kenney (1983) and Dauchot and Dath (1984) used ignition points to fit a kinetic expression. Lyberatos et al. (1984) suggested the use of bifurcation caused by periodic forcing for kinetic parameter fitting. Conrad and Treguer-Seguda(1984) used bifurcation points for parameter estimation of immobilized enzyme systems. The purpose of this study is to examine the advantages and disadvantages of estimating kinetic parameters from a set of experimental ignition and extinction points. 1987 American Chemical Society
Ind. Eng. Chem. Res., Vol. 26, No. 10, 1987 2093 Moreover, we wish to examine the ability of a kinetic model that is based on a fit of bifurcation points to predict the steady-state behavior at operating conditions that are well removed from these points. Development of Steady-State Model We consider a single bimolecular exothermic reaction A + vB products
-
which occurs in a porous catalytic pellet. We assume that the pellet temperature is uniform but different from that of the ambient gas. The corresponding steady-state energy and species balances are
where r(Cs, T,) is the intrinsic reaction rate per unit pellet volume, q(Cs,T,) is the effectiveness factor, and vi is a stoichiometric coefficient. We assume without any loss of generality that A is the limiting reactant and define
(-mkcACbA( To)
(3) hT0 where Tois a reference temperature. Assuming ideal gas behavior and constant total pressure, eq 1 and 2 can be combined to give
P=
4 f f
zA=1+-
-Y) P
VBkcACbA(Tg)(1 - ZA) zB=1-
kcBCbB(Tg)
(44
2
2
N
s = k = l lCaklC(q - @ ) k j ( q - b)lj =l j=1
(9)
where N is the number of observed bifurcation points and q and 4 are the predicted and measured temperature at the bifurcation points, respectively. Subscripts 1 and 2 refer to pellet and gas temperature, respectively, and bkl are elements of the inverted error variance-covariance matrix (Z)-l, where {Zkl! = E ( c k c l ) (10) t k denotes the q - 0 error and E denotes the expected value. A standard, nonlinear, multiresponse optimization program can be used to find poptwhich minimizes S. We used a modified Levenberg-Marquardt algorithm (Froment and Hosten, 1984), which is described in the Appendix. Even when the correct functional form of the rate expression is known from an analysis of the qualitative multiplicity features (Harold et al., 1987a,b), the convergence of the optimization routine depends on the availability of reasonable initial estimates of the unknown kinetic and transport parameters. An initial estimate of some of the parameters may be obtained from T p- Tgvalues along the kinetically controlled (low-temperature) branch. We note that according to (4a)
(4b)
These relations can be substituted into the energy balance to give a single dimensionless steady-state equation
F b , a , P) = 0
Parameter Estimation Scheme A simultaneous solution of ( 5 ) and (8) for any specified p determines y and a or, equivalently, T pand Tg.Clearly, for some p , these equations have no solution; i.e., no bifurcation point exists. The parameter estimation is based on finding p which minimizes the sum of squares
(5)
where p is a vector of kinetic and transport parameters. For example, in the case of a first-order reaction with an Arrhenius temperature dependence,
where
I t is assumed in this derivation that the transport coefficients are independent of temperature and composition. This assumption is relaxed later. The diagram of y versus a is referred to as a bifurcation diagram. When the bifurcation parameter a is changed slowly, the number of solutions of the steady-state equation changes at the bifurcation points, which can be found by a simultaneous solution of eq 5 and
The projection of the bifurcation points in the y direction, obtained by eliminating y between (5) and (8), is called the bifurcation set. A bifurcation m a p is a two-dimensional cross section of the bifurcation set.
where
(-mkcAcbA(Tg) 2 max ( T , - T,) (12) h This adiabatic temperature rise can be determined from experiments at high pellet temperatures. Diagrams of T p - Tgversus CsAcorresponding to several Tp’smay be used to get an initial estimate of the functional dependence of the rate on surface concentration (Le., reaction order) and pellet temperature (activation energy). In addition, the measured AT,, can be used to estimate kcA/h. The mathematical model may predict that a unique solution exists for some set of parameters p for which multiplicity is observed experimentally. A method of overcoming this difficulty is described in the Appendix. ATa, =
Impact of Ambient Temperature Changes The steady-state eq 6 accounts for changes in the total concentration with bulk temperature for experiments conducted under a constant pressure. In general, this expansion effect is not accounted for, and we shall investigate its impact before analyzing the experimental data. To simplify the treatment, we assume that no intraparticle diffusion limitations exist, i.e., 1) = 1. When the impact of the temperature on the total concentration is ignored, the steady-state eq 6 becomes a - y + DUX@ + a - y ) = 0 (13) The highest order singular point of eq 13 is a cusp. Only two possible types of bifurcation diagrams of y versus a exist: a monotonic dependence or one with a counterclockwise hysteresis loop. Moreover, a single a value corresponds to any specified surface temperature y . On
2094 Ind. Eng. Chem. Res., Vol. 26, No. 10, 1987
---
-
P
0
Model I Model II
Ignition
c" Region (a)
!?
1
200
m L a,
Ta
-
0
0
.25
.50
Q
75
c
Da Figure 1. Cross section of the hysteresis and isola varieties of (6) for y = 20 and the Tpversus Tgbifurcation diagram in each of the three regions (insets).
the other hand, inspection of the steady-state eq 6 indicates that two or no a values correspond to any y value. Moreover, (6) has a family of pitchfork singular points at CY = 0.257 y = 0.57 Da = 0 p (14)
-
These points satisfy the conditions F = aF/ay = a2F/ay2 = aF/acu = o
(15)
To divide the global parameter space into regions with different bifurcation diagrams, we follow the scheme used by Balakotaiah and Luss (1982) and construct the hysteresis (H) and isola (I) varieties, which satisfy the conditions H:
F = dF/dy = d2F/dy2 = 0
(16)
I:
F = aF/ay = aF/acu= o
(17)
Using (6) and (16), we find that the hysteresis variety is defined by
p = - 4cu3
y - 4ff
with cu E (0,0.25y). Using (6) and (171, we find that the isola variety is given by D a =yT- e4cux p [
Y (12cu -2 4
]
P=-
a2y
y
-
4ff
(19)
Figure 1 describes these varieties in the (Da, 6 ) plane for y = 20. The varieties divide the feasible region Da L 0 into three different regions with three qualitatively different bifurcation diagrams of T , versus Tg(or y versus a ) , shown by the insets. These diagrams are rather different from those predicted by the commonly used steady-state eq 13. First, at sufficiently low gas temperatures, a branch exists along which T, is a monotonic decreasing function of Tg.This feature cannot be predicted by (13) for any exothermic reaction and has not been observed in previous experimental studies. One reason may be that the limited range of gas temperatures may prevent one from observing the negatively sloped high pellet temperature branch. For example, in region b, one finds only the common hysteresis if the lowest gas temperature is not sufficiently small. Second, a surprising finding is the existence of isolated branches such as those shown in case
(0
m
(3
4
5
6
XbE
(%v. )
7
Figure 2. Experimental and model predictions of the boundaries of the multiplicity region. Schematic bifurcation diagrams of T, versus Tg(xbE)are shown in circles (squares).
c in Figure 1. In such a case, no extinction of the upper state is obtained following a continuous slow decrease of the feed temperature. In practice, the feed temperature is larger than some minimum so that lack of an extinction is not sufficient to discriminate between cases b and c. The hysteresis and isola varieties shown in Figure 1 intersect at a pitchfork point located on the boundary of the feasible region (Da = 0, p = a). Another type of bifurcation diagram exists for Da < 0. This example indicates that the change in total concentration with ambient temperature at constant pressure may have an important impact on the multiplicity features, and it is essential to account for it in any fit of experimental data.
Analysis of Experimental Data: Oxidation of Ethane The proposed procedure is illustrated by applying it to experimental data obtained during the oxidation of ethane on a single 4-mm-diameter 0.3% wt Pt on alumina spherical pellet. Most of the platinum was impregnated within an exterior shell of about 100-pm thickness. The pellet was immersed in a stream of ethane, oxygen, and nitrogen. A detailed description of the experimental setup and procedure was presented by Harold and Luss (1985). The experiments were conducted at 1-atm total pressure over a gas temperature range of 20-300 "C and ethane concentrations in the range 3.63-5.81% vol. The linear gas velocity was about 1.55 cm/s. A unique stable state existed for sufficiently high or low gas temperatures, and no more than two stable states were found under any operating condition. Two types of bifurcation diagrams of Tpversus Tgwere observed, a single-valued monotonically increasing curve and a counterclockwise hysteresis loop. The ignition and extinction gas temperatures decreased with increasing reactant concentrations and the region of two stable states was bounded within a cusp (Figure 2). No partial oxidation products were found during the oxidation of low molecular weight alkanes, olefins, and alcohols on noble metal catalysts in the studies of Hiam et al. (1968), Moro-oka and Ozaki (1966), Moro-oka et al. (1967), Schwartz et al. (1971), and Yao (1980). The only exception was the oxidation of 2-propanol on Pt (Schwartz et al., 1971) in which acetone was obtained under certain
Ind. Eng. Chem. Res., Vol. 26, No. 10, 1987 2095 at each T pwas fit to the empirical relation -
A
280
23
0
32C
48
x,E
(%
V.
T p- T g =
The scatter of the data is due to the approximate analysis and experimental error. The temperature dependence of b satisfies an Arrhenius relation with an apparent activation energy of 61.4 kJ/mol. These estimates of n, E , and kc/h are very useful initial guesses in the parameter fit algorithm. We assume that under the experimental conditions, heat transfer by radiation can be neglected. Following the Ranz and Marshall (1952) correlation and the Chilton and Colburn (1934) analogy for heat and mass transfer, we assume that
CpH,)
Figure 3. Dependence of temperature rise on estimated ethane surface concentration a t several fixed surface temperatures. Solid line represents least-squares fit.
operating conditions. Thus, we assume that the oxidation can be described by a single overall reaction C2H6 3.502 2C02 + 3H20 -AH = 1429 kJ/mol (20) To estimate the impact of intraparticle diffusion resistance, we applied the Weisz-Prater criterion (1954) which for a pellet with an active exterior shell in R1 < R < Rp is defined as
+
-
The calculations showed that this modulus was much smaller than unity, implying that the effectiveness factor is very close to unity. Estimates of the maximum intraparticle gradients showed that they were always smaller than 10% of the external temperature gradients. We therefore neglected the intraparticle gradients. Previous studies of ethane oxidation on a Pt wire (Hiam et al., 1968) or Pt/A1203catalysts (Yao, 1980) reveal that the isothermal rate increases monotonically with surface ethane concentration. The cusp-shaped bifurcation map (Figure 2) is indicative of a rate which increases monotonically with an increase in the concentration of the limiting reactant (Harold et al., 1987a). Yao (1980) claims that the rate has a negative order (-3) with respect to oxygen. Thus, despite the expected complex reaction mechanism, we assume that the overall rate can be described by the empirical expression (22)
where subscripts E and 0 denote ethane and oxygen, respectively. The adiabatic temperature rise, AT,,, was determined from experimental data at high pellet temperatures by using (12). The value of the adiabatic temperature rise per 1% ethane in the gas phase, ATad/lOOXbE,was found to be 50 Oc/1% vol ethane, where XbE is the bulk ethane mole fraction. By use of this value, (12) predicts that kcE/h is equal to 155 K cm3/J for = 0.5(T + T,)= 679 K. An estimate of the functional form of t i e isothermal rate was obtained by plotting the temperature rise at several constant pellet temperatures versus ethane surface mole fraction ( x , ~under ) constant total pressure. (The surface mole fraction is computed by the expression x S E = zExbETp/Tg, where z E is determined by (12).) Figure 3 shows that the isothermal rate increases monotonically with X,Efor pellet temperatures in the range 260-320 "C for which the reaction is kinetically controlled. The rate
h = h(To)(T'/T0)0'25
(244
kc = k,(To)(T/T0)1.25
(24b)
T' = wTg + (1- w)TP
(25)
where is a weighted average of the gas and pellet temperatures. We discuss later the proper choice of the weight factor, w. Substitution of the rate expression given by (22), the ideal gas law for Chi, and the transport coefficient temperature dependences (eq 24a and 24b) into eq 1and 2 and defining
D=
dC,, T,)= ~ ( T , ) C , E " C , O ~
(23)
~ X , E ~
k (To)[ c b ~ ( T o ) [&(To) 1
vp
kcE( T0)Sx
where
The experiments were carried out using a large excess of nitrogen so that the ethane-oxygen mass-transfer coefficient ratio was computed by the expression kcE/kcO De,EN2/De,0Np (30) and was estimated to be about 0.77. Some comments on the sensitivity of the optimum to this diffusivity ratio are provided later. The steady-state eq 27 contains six unknown parameters (n, m, E , k,E(TO)/h(TO),w , and ko/ k c E ( TO)) The values of (y, a) at the bifurcation points for any p were computed by a simultaneous solution of (27) and
YzEzO
Y2
[ w a + (1- w ) y ] -
(a - Y ) ( l - w ) w a (1- w ) y
+
.
1.
2096 Ind. Eng. Chem. Res., Vol. 26, No. 10, 1987 Table I. Optimal Parameter Values for Models I and I1 and the Corresponding 95% Confidence Limits and Sum of Squares h / MTo), ( ~ m ~ ) " +I~(g-mol -l/ ~,E(TO) / h(To), model n m u: E E)"-'(g-mol Oz)"Ycmpellet)] cm3 K / J S I 4.55 f 1.49 0 0.7 130.6 i 32.1 6.56 x 1031 f 1.74 x 1033 375.1 f 11.0 552 I1 3 -3 0.4 101.5 i 1.09 8869 f 230 265.6 i 4.3 709 - -
-\--\
600
h
300 Model I
a,
Made1 I1
200
3
5
4 XbE
6
(%v.)
- -- -
i
1 ooo
100
200
300
G a s T e m p e r a t u r e , Tg ("C)
Figure 4. Comparison of measured pellet temperature a t ignition and extinction with the predictions of the two models.
Figure 5. Predicted bifurcation diagrams of Tp versus T8 for a mixture with 6% vol ethane.
Simulations were carried out for two cases. In model I, it was assumed that due to the large excess of oxygen ( C b O / C b E exceeded always 3.5, the stoichiometric ratio), the rate is independent of the oxygen concentration; Le., m = 0. In model 11, the values of n and m were assigned integer values to simplify the parameter estimation, and the sum of squares was used to find the combination of n and m which gave the smallest value. For an assumed n and m, the value of the weight factor (w)was assigned different integer multiples of 0.1 between 0 and 1. The w value which gave the best fit (i.e., lowest S ) for a particular n and m combination was taken as the correct estimate. The parameter estimation scheme requires knowledge of the measurement errors. Based on our experience, we estimate that the pellet and gas temperature variances are 9 and 4 K2, respectively, and that the pellet and gas temperatures are uncorrelated. Table I reports the estimated parameters and the corresponding 95 % confidence limits and sum of squares for both models. The first model predicted a very high reaction order for ethane (4.55). The rather poor confidence limits for most of the parameters are indicative of the strong correlation between these values. The multidimensional surface of sum of squares was very flat next to the minimum, and convergence was slow. The sum of squares for the second model was slightly larger than that for the first model, as it contains only three parameters. However, the confidence limits of the parameters are significantly reduced. Moreover, the predicted values of n = 3, m = -3, and E = 101.3 kJ/mol agree closely with those reported by Yao (1980) ( n = 2, m = -3, E = 83.7 kJ/mol). Equation 24 predicts that the estimated k,E(To)/h(TJ of 265.6 cm3 K / J at To = 1000 K is equal to 180.6 at 679 K, which is close to 155 K predicted by the adiabatic temperature rise estimate discussed earlier. Figure 4 compares the measured pellet temperatures at the bifurcation points with values predicted by the two models. This fit is inferior to that of the gas temperatures (Figure 2). The main differences between the predicted values and the data are along the branch of extinction points. Model I1 fits the diagrams better than model I. Both models I and I1 predict, as does the model for a first-order reaction discussed earlier, that an isola point
exists when the gas temperature is the bifurcation variable. Schematic Tp versus Tgbifurcation diagrams are shown in circles in Figure 2, while those of Tpversus XbE are shown in squares. Clearly, four types of Tpversus Tgdiagrams exist. Crossing the isola point causes a transition from type c to d. In practice, parts of these diagrams may not be observable, as they correspond to unfeasible operating conditions. The dashed vertical lines denote the feasible boundary for T > 0 " C and XbE C 6.0% vol (arbitrary). Thus, the only dfference between cases b and c is that the first has two bifurcation points in the feasible region, while c has three. The model predicts that the extinction loci approaches asymptotically Tg= 0 K and Tp m as XbE 0. This is an erroneous prediction, as the ideal gas law, and consequently the model, is not valid at very low gas temperatures. Two types of Tp vs. Tgdiagrams were observed experimentally (shown as cases a and b in Figure 2). The simulations indicate that for XbE larger than 5.85% vol, a bifurcation diagram of type d exists. Figure 5 shows such a case for XbE = 6% vol. We were not aware of this point a t the time and carried out experiments only for XbE 5 5.81% vol. This upper bound was chosen to keep the temperature well below 600 "C to minimize sintering. Unfortunately, the original pellet was destroyed by the time the simulations revealed this interesting feature. Experiments with a similar pellet and XbE = 8% vol gave a bifurcation diagram similar to that shown in Figure 5, i.e., no extinction point, but failed to reveal a minimum in the pellet temperature along the ignited branch. The inability to observe such a minimum was caused by the experimental constraint of Tg> 0 "C. Thus, we cannot conclude if an isolated branch existed. The simulations show that the value of the weight factor (w), which is needed to account for the temperature effect on the transport coefficients, has a strong impact on the quality of fit. For example, by use of w = 1, the predicted s;2" at extinction are much larger than those observed for all XbE 6% vol. In contrast, for w = 0, both models predict that an isola point exists for 3tbE well within the experimental region (15.81% vol), even though it was not observed experimentally. Intermediate values of w improve significantly the fit. In addition, the optimum was
-
-
Ind. Eng. Chem. Res., Vol. 26, No. 10, 1987 2097 250
300
350 b . 4.56%
. O
350
350
250
250
1
2
E
300
200 200
250
3000
100
Gas Temperature,T,
200
100 300
('GI
Figure 6. Comparison of four experimental bifurcation diagrams of Tp versus Tgwith predictions of the two models.
rather sensitive to the assumed diffusivity (mass-transfer coefficient) ratio (eq 30). When the value of 0.77 was varied within & l o % , it was found that the fit may be improved by adjustment of this value. However, the variation is limited physically, so we retained 0.77. An important objective of this study was to examine the quality of the fit of the entire bifurcation diagrams by a kinetic model, the parameters of which are based only on a fit of the bifurcation points. Figure 6 compares the experimental bifurcation diagrams of T pversus Tgfor four levels of XbE in (3.63, 5.81% vol) with the predictions of the two models. The diagrams indicate that while the predictions of the bifurcation points are adequate, those of the bifurcation diagrams are much less satisfactory. A t low pellet temperatures, both models underpredict the measured values, while along the ignited branch the temperatures are overpredicted. In general, the predictions of the bifurcation diagrams by the second model are more accurate than those of the first model. The overestimate of the pellet temperature of the ignited states is probably caused by radiation heat losses, which are not accounted for in the model. Obviously, this error increases with increasing pellet temperatures.
Concluding Remarks This work is the first attempt to use an observed bifurcation map to fit the parameters of a kinetic model of a catalytic reaction. Several important conclusions can be drawn from this study. The simulations show that it is essential to have a reasonable initial guess of the parameters so that the model predicts multiplicity at least for a large fraction of the cases in which it was observed. The main difficulty is guessing the kinetic and transport parameters, and this may be done by a preliminary analysis of the temperature rise along the low-temperature branch and high-temperature branch, respectively. Our results indicate that secondary effects, which are often ignored in the literature, can have an important impact on the bifurcation map features and have to be accounted for in any fit of experimental data. The strong sensitivity of the fit to the value of the weight factor ( w ) and the ratio of oxygen to ethane mass-transfer coefficients is rather surprising. The uncertainty involved in the values of these two parameters affects adversely the quality of the data fit. The existence of an isola point due to the impact of changes in the bulk temperature on the bulk
concentration is a novel finding, which has not been reported so far. The results suggest that the quality of the fit can be improved if one fits both the bifurcation map as well as several bifurcation diagrams, i.e., nonsingular points. Specifically, one finds the set of parameters which minimizes the sum of squares S' = + tS(T,, T J (32)
s
where t is a weight factor and S(T,, T,)is the sum of squares of the difference between the measured and predicted gas and pellet temperatures for M nonsingular points on the bifurcation diagrams. Clearly, when t = 0, S' = S. It is desirable to start the fit with a relatively large value of t to minimize the impact of a poor initial guess on the value of S. The value of t can then be decreased as the estimate of the parameters improves. We expect this fitting procedure to give a better fit than the one which is based only on the bifurcation points. The predictions of both models break down for sufficiently low Tgat which the ideal gas law is no longer valid. This deficiency may be overcome by use of another equation of state, but this will be of little value, as one is not expected to conduct experiments in this range of very low gas temperatures.
Acknowledgment We are grateful to the National Science Foundation, the Welch Foundation, and the Shell Companies Foundation for support of this research.
Nomenclature Cb = bulk concentration C, = surface concentration D = Damkohler number defined in (26) Da = Damkohler number defined in (7) De = effective diffusivity E = activation energy F = function defining steady state h = heat-transfer coefficient AH = heat of reaction I = identity matrix k = rate constant k , = mass-transfer coefficient lzo = preexponential factor m = reaction order M = number of nonsingular points M = sensitivity matrix defined in (A4) n = reaction order N = number of experiments (Le., bifurcation points) Nd= number of concentrations for which multiplicity is observed Np = number of estimated parameters p = vector of parameters q (4) = experiment (calculated) temperature r = reaction rate c = vector of residuals R = radial coordinate R, = gas constant R, = pellet radius R1 = distance from pellet center to active shell S = sum of squares S , = external pellet surface area S' = modified sum of squares t = weight factor for modified sum of squares T = temperature T = gas temperature = reference temperature T = pellet temperature r"= average film temperature ATad = adiabatic temperature rise V , = pellet volume ~
2098 Ind. Eng. Chem. Res., Vol. 26, No. 10, 1987
IrUh= (qlu - hu)p=pz1 = 1, 2
w = weight factor for film temperature xb = bulk mole fraction n, = surface mole fraction X = defined in (7) y = dimensionless pellet temperature z = dimensionless surface concentration
where the index u refers to a bifurcation point and u to a parameter. The elements of M are determined by an analytical or numerical differentiation of (6) and (8) with respect to p u and solving (Mu,,lkifor k = 1 and 2. The fitting procedure encounters difficulties when for a given set of parameters, multiplicity is not predicted for Nd experiments for which it was observed. In these cases, some of the elements in both Mkiand rliare not defined. We overcame this problem by replacing the nonexisting elements of Mkiby those computed in the previous iteration and the Nd missing elements of rliby the average residual
effectiveness factor
y = dimensionless activation energy A = adjustable stepping parameter defined in (A21 Y = stoichiometric coefficient 4 = Weisz-Prater modulus u
= element in inverted variance-covariance matrix
Z = variance-covariance matrix Subscripts ad = adiabatic A = chemical species b = bulk
c = for mass-transfer coefficient e = effective E = ethane 6 = pas i = zth element, species, etc. j = jth element, species, etc. K = kth element, species, etc. 1 = lth element, species, etc. N2 = nitrogen 0 = molecular oxygen p = pellet s = surface u = uth element, species, etc. u = uth element, species, etc. x = external 0 = reference
This scheme predicted the parameters well for cases in which Ndwas a small number. A more detailed description of the scheme was presented by Harold (1985). Registry No. C2Hs, 74-84-0.
Literature Cited
Superscripts i = ith iteration T = transpose
Appendix: Numerical Scheme To Find Minimum of s One starts by assuming a set of parameters p and computes the corresponding bifurcation gas and pellet temperatures for each of the N experiments, by solving simultaneously (6) and (8). The value of S is then computed by (10). If a specified convergence criterion is not satisfied, the values of the parameters are updated by the formula = pi 2
Api =
2
+ A@' 2
(AI)
2
C C(aklMk?Mli+ XI)-'C CuklMk?rli
k=l1=1
(A3)
Mki is the N X Np(number of parameters in p ) sensitivity matrix, the elements of which are IMuulki = ( @ k u / d ~ u ) p = p r k = 1, 2 (A4) u = 1,...,Np u = 1,...,N
Greek Symbols a = dimensionless gas temperature p = dimensionless temperature rise t = error 9 =
u = 1,-.,N
k=l1=1
(A2)
where rliis the N-dimensional residual vector, the elements of which are
Balakotaiah, V.; Luss, D. Chem. Eng. Sci. 1982,37,1611. Beskov, V. S.;Karavaev, M. M.; Zharov, D. V.; Arutyanyan, V. A. React. Kinet. Catal. Lett. 1976,4 , 351. Chilton, T. H.; Colburn, A. P. Ind. Eng. Chem. 1934,26,1183. Conrad, F.;Treguer-Seguda, V. Chem. Eng. Sci. 1984,39,705. Dauchot, J.-P.; Dath, J.-P. J . CataE. 1984,86,373. Froment, G.F.;Hosten, L. In Catalysis Science and Technology; Anderson, J. R., Boudart, M., Eds.; Springer-Verlag: New York, 1984; Vol. 6,p 97. Harold, M. P. Ph.D. Dissertation, University of Houston, 1985. Harold, M. P.;Luss, D. Chem. Eng. Sci. 1985,40,39. Harold, M. P.; Sheintuch, M.; Luss, D. Znd. Eng. Chem. Res. 1987a, 26, 786. Harold, M. P.; Sheintuch, M.; Luss,D. Znd. Eng. Chem. Res. 1987b, 26, 794. Herskowitz, M.; Kenney, C. N. Can. J . Chem. Eng. 1983,61, 194. Hiam, L.; Wise, H.; Chaikin, S. J . Catal. 1968,9,272. Jiracek, F.;Havlicek, M.; Horak, J. Collect. Czech. Chem. Commun. 1971,36,64. Lyberatos, G.; Kuszta, B.; Bailey, J. E. Chem. Eng. Sci. 1984,39,739. Moro-oka, Y.; Ozaki, A. J . Catal. 1966,5, 116. Moro-oka, Y.; Morikawa, Y.; Ozaki, A. J. Catal. 1967,7, 23. Rader, C. G.;Weller, S. W. AIChE J. 1974,20, 515. Ranz, W. E.; Marshall, W. R., Jr. Chem. Eng. Prog. 1952,48, 141. Schwartz, A,; Holbrook, L.;Wise, H. J . Catal. 1971,21, 199. Weisz, P. B.; Prater, C. D. Adv. Catal. 1954,6,143. Yao, Y.-F. Ind. Eng. Chem. Prod. Res. Dev. 1980,19,293. Received for review March 17, 1986 Accepted June 1, 1987
