433
Ind. Eng. Chem. Fundam. 1983, 22, 433-436
n = number of group increments (use in Table I) ASi = group increment for eq 2; cal mol-’ K-’ Asvb = entropy of vaporization at the normal boiling point,
cal mol-’ K-’ Tb = normal boiling point, K Subscripts cal = calculated value lit. = literature value ASvb = Z A S j =
E(-CH,)
+ E(-CH2-) +
Literature Cited
nonring incr
E(=CH-)
+ Z(=C) + E(doub1e bond number between carbon atoms) ring incr
+ + + 0.3437 X 7)+ (-9.3621 X 3) + (-27.1977 + 9.1993 X 5)= 21.6945 cal mol-’ K‘’
= (9.5937) (0.2320) (18.7504
error: 0.3%. Nomenclature E = error, % AHvb= latent heat of vaporization at the normal boiling point, cal mol-’ Ncompd = number of compounds
Chen. N. H. J. Chem. Eng. Data 1085, 10, 207. Giacalone, A. Gazz. Chim. Ital. 1051, 8 1 , 180. Hoshlno, D.; Nagahama, K.; Hlrata. M. J. Jpn. Pet. Inst. 1070, 22, 32. Shinoda, K. “Principles of Solution and Solublltty”; Marcel Dekker: New Ywk, 1978; p 8. McCurdy, K. G.; Laldler, K. J. Can. J . Chem. 1063, 4 1 , 1867. Miller, D. I d . Eng. Chem. 1084, 56(3), 46. Procoplo, J. M.; Su, 0. J. Chem. Eng. 1068, 74(12). 101. ReM, R. C.; Prausnk. J. M.; Sherwood, T. K. “The Propertiis of Gases and Liquids”, 3rd ed.; McGraw-Hlll: New York, 1977; Chapter 2 and 6. Rledel, L.; Chem. Ing. Tech. 1054, 26, 679. Rosenbrock, H. H. Compuf. J . 1080. 3 , 175. Wllhok, R. C.; Zwollnskl, B. J. “Handbook of Vapor Pressure and Heats of Vaporization of Hydrocarbons and Related Compounds”, American Petroleum Institute Research Project 44, Texas A and M University, College Station, TX 1971; Sectlon 111.
Received for review August 17, 1982 Revised manuscript received May 23, 1983 Accepted July 8,1983
Aspects of Autocatalytic Reaction Kinetics Michael Frenkiach’ and David Clary Department of Chemlcal Engineering, Louisiana State Unlvers& Baton Rouge, Louisiana 70803
The kinetics of the autocatalytic system (R1)-(Rz) is investigated. A notable feature of this system appears when the dependency of the concentration of product C on the ratio r = k,/kzCAO at a given time is considered. Analysis shows that the effect of variation in r on the concentration of C depends upon the manner in which r is altered. A maximum in Cc vs. r is observed when the numerator of r is varied while the denominator is held constant; however, no maximum is observed when only the denominator is varied. Another noteworthy feature of this system is the existence of “frozen” composition at infinite reaction time.
Introduction The authors of most discussions of autocatalytic reaction kinetics (Benson, 1960; Froment and Bischoff, 1979; Hill, 1977; Houser, 1969; Levenspiel, 1972) usually concentrate on the time development of the system. This is probably due to the fact that the basic characteristic of an autocatalytic reaction is the presence of a maximum in reaction rate as a function of time or reactant concentration. However, some other fundamental aspects of the autocatalytic reaction appear to be of theoretical interest and practical importance. For example, the analysis of these aspects for the reaction sequence kl
A-B
(RU
reactant A to final product C. This or a similar sequence is a typical mechanism for a variety of natural processes for which reaction R1 (presumably a decomposition reaction) has a much greater activation energy than does reaction R2. This paper presents a detailed analysis of the reaction sequence introduced above.
Mathematical Development The kinetics of the system (Rl)-(R2) is described by the following set of differential equations dCA/dt = -k,CA - ~ & A C B (la) dC,/dt = klCA - k2CACB dC,/dt = ~ ~ C A C B
k2
A+B-C
(R2) has been successfully applied to the global kinetics of soot formation during the pyrolysis of aromatic hydrocarbons (Frenklach et al., 1983). The above sequence (Rl)-(R2) is an autocatalytic process with respect to the intermediate product B, while the overall reaction is the conversion of
with the initial conditions CAlt=O =
c&
CSl,=O = 0 C&o = 0
0196-4313/83/1022-0433$01.50/00 1983 American Chemical Society
(1b) (IC)
434
Ind. Eng. Chem. Fundam., Vol. 22, No. 4, 1983
Expressed in dimensionless variables, eq la, lb, and IC become (24 da/dt‘= -a - ab/r db/dt‘= a - ab/r
(2b)
dc/dt‘= ab/r
(2c)
with the initial conditions alt,=o= 1
blt,=o = 0 CJtkO =
0
where a = CA/C,, b = C,/C,, c = Cc/C,, t’= kit, and r = kl/lz2CA,. It is not obvious that a complete analytical time-dependent solution of the set of eq 2 is possible. However, using the elegant substitution (Holland and Anthony, 1979) d7 = k,a dt
(3)
1
log
r
Figure 1. Dimensionless concentration of product C vs. log r, where the ratio r = kl/kZCAO is varied by changing kl while keeping kZcA0 constant. The characteristic time t* in this case is given as t* =
the nonlinear differential eq 2 are transformed into linear equations (44 d a / d r = -1 - b/r db/dr = 1 - b/r
(4b)
dc/dr = b/r
(44
with initial conditions =1 bl,=o = 0 CI,=O
=0
The solution of the set of linear differential eq 4 is readily obtained a = 1 - 27
+ r(1 -
= 1 - 27
+b
(5a) 1
b = r(l c = 7 - r(1 - e - r / r ) =
(5b) 7
-
b
(54
The independent variable of the eq 4 is “transformed time” T rather than real reaction time t. Integrating eq 3, one obtains a relation between these two times
log
r
Figure 2. Dimensionless concentration of product C vs. log r, where the ratio r = kl/k2CAois varied by changing kzCAowhile keeping kl constant. The characteristic time t* in this case is given as t* = l / k l .
constant, the reaction time t may be expressed in terms of the characteristic time tl* = l/k,CAo; that is
t = constant
X
tl*
(8)
Substituting (8) into (7), one obtains Substitution for a(T) from (5a) yields
For a given value of r, the integral on the right-hand side of eq 7 can be solved numerically. Thus, the system composition for the chosen value of r may be found as a function of the reaction time t.
Discussion The time development of the reaction system (Rl)-(R2) (i.e., the time-dependency of the species concentrations) exhibits behavior similar to that analyzed in detail elsewhere (see, e.g., Levenspiel, 1972). An interesting feature appears when one considers the dependency of the concentration of product C on the ratio r at a fixed reaction time. The dependence of Cc on r is changed according to whether r is varied by changing ita numerator kl or its denominator kzCAo. If k , is varied while kzCAo is held
For a given value of r (i.e., for a given value of kl in this case), the numerical solution of eq 9 determines T and hence the concentrations of A, B, and C at a given time t. The computational results indicate that Cc/CA0passes through a well-defined maximum when plotted vs. r as k , is varied (Figure 1). If, however, kl is held constant and the ratio r is varied by changing k2CA0,the appropriate time scale is tz* = l / k l and the reaction time may be expressed as t = constant’ X tz* (10) Substituting (10) into (71, one obtains constant’ =
dr 1 - 27
+ r(1 - e-r/r)
(11)
:/
0
...
s-
0
20 -
0
4
-
0
Ind. Eng. Chem. Fundam., Vol. 22, No. 4, 1983
I
’9-
-
435
I
Y
-0.
\
E-
W
-
N-
9 p3.0
-2.0
-1.0
1.0
0.0
log
2 ..‘o0
0
r
Figure 4. Dimensionless concentrations of products B and C at infinite reaction time vs. log r.
O 0
+.UO
1:s:
1:BO
2’.00
TEMPERRTURE
2’.20
2’.UO
I
2.60
[KI / l o 3
Figure 3. Soot yield versus temperature at different reaction times: ( 0 ) 0.5 ms; (A) 1.0 ms; (0) 1.5 ms. The solid lines are cubic spline approximations to the data. The results were obtained in the pyrolysis of a 1.75% toluene-argon mixture (total carbon atoms concentration -1.7 X 10’’ atoms/cms) behind reflected shock waves by monitoring the extinction of a 632.8-nm H e N e laser beam (Frenklach et al., 1983).
The numerical solution for this case shows that there is no maximum in product concentration vs. the ratio r (Figure 2). Mathematically, the difference between the two cases is apparent when expressions (9) and (11)are analyzed. Phenomenologically, this difference in behavior may be explained in the following manner. When reaction R1 is slow, reaction R2 is also slow due to the low concentration of intermediate B. Increasing the rate constant kl will enhance the production rate of B, which, in turn, increases the product of concentrations, CACB, thus increasing the amount of C formed. However, a further increase in kl will begin to influence the product term, CACB, by lowering the concentration of A; this will retard the formation of C. Therefore, when compared at the same reaction time with k2 held constant, the concentration of product C must pass through a maximum as kl is increased. On the other hand, if kl is held constant, increasing k2 will enhance the rate of interaction between A and B in a monotonic manner and therefore no maximum in Cc vs. r will be observed. The reaction sequence discussed in this paper may be useful as a global model for some natural phenomena whose mechanisms have not been well-established. For example, in a soot formation model proposed by Frenklach et al., reaction R1 represents fragmentation of an intact aromatic ring which is presumably a high activation energy process. Reaction R2 is assumed to model a free-radical polymerization and therefore is expected to have a relatively low activation energy. An increase in temperature would increase both kl and k 2 but, due to the relative magnitudes of the activation energies, the ratio r would increase with temperature. Figure 3 represents the typical experimental results of soot yields obtained at various temperatures (Frenklach et al., 1983). The similarity of Figure 3 and Figure 1 is apparent. Another distinctive feature of the autocatalytic sequence
(Rl)-(R2) is the existence of “frozen” composition at infinite reaction time. Figure 4 depicts the dependency of the concentrations CBmand Cc., on the ratio r (it should be noted that these concentrations are functions only of the ratio r and not of how r is varied). This feature of kinetic behavior may be significant. Typically, when a researcher studying a poorly understood reaction system in a batch reactor observes that the product concentrations have reached a time plateau, he assumes that equilibrium has been attained. The above analysis indicates that an alternative explanation for experimentally observed time-plateaus in concentrations attained in a batch reactor may be the assumption of “frozen- composition if the process may be globally modeled as (Rl)-(R2) or a similar reaction sequence. Conclusions The analysis presented in this paper shows that the autocatalytic sequence (Rl)-(R2) exhibits some peculiarities. The concentration of product C, when plotted vs. the ratio r at a given reaction time, passes through a well-defined maximum if the rate constant k , is varied while k2CA0is held constant. However, no maximum exists when k2CAO is changed while kl is held constant. The possible utilization of these results follows from the fact that the sequence (Rl)-(R2) is usually a global description of real phenomena and therefore the parameters kl and k2 may depend on a variety of factors. Another distinctive feature of the reaction sequence (Rl)-(R2) is the existence of “frozen” system composition at infinite reaction time. This result indicates that the experimental observation of concentration time-plateaus in a batch reactor may not necessarily mean that equilibrium has been achieved but may instead suggest the possibility of a “frozen- composition. Acknowledgment D.C. wishes to acknowledge support from the National Science Foundation’s Graduate Fellowship program. The work was supported by the Office of Fossil Energy, U.S. Department of Energy, under the auspices of Grant Number DE-FG22-80PC30247. Publication costs assisted by NASA Contract No. NAS3-23542. Nomenclature t = reaction time t’= kit t* = characteristic time C,,[A] = concentration of reactant A
Ind. Eng. Chem. Fundam. 1983,22, 436-445
438
CB,[B] = concentration of intermediate product B Cc,[C] = concentration of final product C CAo,[Al0 = initial concentration of reactant A CB,,[B]- = concentration of product B at infinite reaction time Cc,,[C], = concentration of product C at infiiite reaction time a = CA/CAO =
CB/CAO
c = cC/cAO
= kl/k2CA0 k1 = rate constant for reaction R1 k2 = rate constant for reaction R2 E ( m ) = -Im[(m2 - l)/(rf + 2)] m = complex refractive index of soot particles T = transformed time, defined by eq 3
Literature Cited Benson. S. W. “The Foundations of Chemical Kinetics;” McGraw-Hill: New YO&, 1960;pp 45-46. Frenklach, M.; Taki, S.; Matula, R. A. Combust. Flame 1983, 49, 37. Froment, G. F.; Bischoff, K. B. “Chemical Reactor Analysis and Design;” Wiiey: New York, 1979;pp 13-15. Hili, C. G., Jr. “An Introduction to Chemical Engineering Kinetics and Reactor Design;” Wiiey: New York, 1977: pp 337-342. Holland, C. D.; Anthony, R. G. “Fundamentals of Chemical Reaction Engineering;” Prentice-Hail: Englewood Cliffs, NJ, 1979;Chapter 10. Houser, T. J. J. Chem. Wys. 1969,50. 3962. Levenspiel, 0. “Chemical Reaction Engineering,” 2nd ed.; Wiiey: New York, 1972;pp 150-157.
Receiued for reuiew September 13, 1982 Accepted July 18, 1983
Improvement of Gauss-Newton Method for Parameter Estimation through the Use of Information Index Nlcolas Kalogerakis and Rein Luus’ Department of Chenflcal Englneerlng and Applled Chemlstty, Unlversky of Toronto, Toronto, Ontario M5S 1A4, Canada
I n parameter estimation for systems described by ordinary differential equations, the region of convergence for the Gauss-Newton method can be substantially enlarged through the use of an “information index” and an optimal step-size policy. The information index provides a measure of the available sensitivity information as a function of time, thereby locating the most appropriate section of data to be used for parameter estimation. The use of the chosen section of data significantly improves the conditioning of the linear algebraic equations that are solved at the end of each iteration in the Gauss-Newton method, and the region of convergence is thus expanded. I f observations are unavailable in the most appropriate section of data, artificial data can be generated by data smoothing and interpolation. The effectiveness of the approach is illustrated by four examples deplctlng typical chemical engineering systems. The proposed procedure is especially attractive for systems described by stiff differential equations.
Introduction Frequently the structure of the mathematical model of a physical system is available from fundamental principles. It is then required to use experimental data to obtain the best estimate for the parameters in the model in the sense of minimizing the difference between the model predictions and the actual observations taken from the physical system. Quasilinearization and Gauss-Newton are the most commonly used parameter estimation procedures because of their fast quadratic convergence near the optimum (Bellman and Kalaba, 1965; Lee, 1968; Bard, 1970,1974). However, a major problem is the small region of convergence. Unless the initial guess of the unknown parameters is in the immediate vicinity of the optimum, divergence occurs. This problem has received much attention recently. Ramaker et al. (1970) suggested the incorporation of Marquardt’s modification to expand the region of convergence. Donnelly and Quon (1970) proposed a procedure whereby the given data are perturbed if divergence occurs and a series of problems are solved until the model trajectories Ymatch”the original data. Nieman and Fisher (1972) incorporated linear programming into the parameter estimation procedure and suggested the solution of a series of constrained parameter estimation problems where the search is restricted in a small parameter space around the chosen initial guess. Recently, Wang and Luus (1980) showed that the use of shorter data length can enlarge substantially the region of convergence. To further in-
crease the region of convergence, Kalogerakis and Luus (1982) proposed a two-step procedure whereby direct search optimization is used initially to reach the vicinity of the optimum, followed by the Gauss-Newton method. Whenever an initial guess is outside the region of convergence the source of difficulties may be traced to the following. First, if the estimate of the parameter vector is far from the optimum, the differential equations may become numerically unstable. Second, the Gauss-Newton method may cause excessive overstepping; and third, it may predict a totally wrong direction in which to change the parameter vector. In this paper the importance of using a robust integration routine is demonstrated and an optimal step-size policy is proposed to cope with the first two problems. It is further shown that lack of sufficient sensitivity information cause8 the Gauss-Newton method to predict a wrong direction and the problem can be overcome by incorporating such information into the parameter estimation procedure. Problem Statement Let us consider a dynamic system described by the ordinary differential equation
x(to) = xo (2) where x is the n-dimensional state vector, xois the given 0 1983 American Chemical Society
