Canjar. L. N., Manning, F. S.,"Thermodynamic Properties and Reduced Correlations for Gases", Gulf Publishing Co., 1967. Chandron, J., Asselineau, L., Rennon, H., Chem. Eng. Sci., 23, 839 (1973). Chang, S. D.. Lu, B. C.-Y., Can. J. Chem. Eng., 48, 261 (1970). Chueh, P. L., Prausnitz, J. M., lnd. Eng. Chem., Fundam. 6, 492 (1967). Das, T. R., Eubank, P. T., Adv. Cryog. fng., 18, 208 (1973). De Mateo, A., Kurata, F., lnd. Eng. Chem., Process Des. Dev., 14, 137 (1975). Din, F. "Thermodynamic Function of Gases", Vol. 3, Butterworths. London, 1961. Elliot, D. G., Chen, R. J. J., Chappelear, P. S.,Kobayashi, R., J. Chem. Eng. Data, 19, 71 (1974). Eubank, P. T., Adv. Cryog. Eng. 17, 270 (1972). Goodwin, R. D., Prydz. R., NBS J. Res., 75A, 81 (1972). Guggenheim, E. A., J. Chem. Phys., 13, 253 (1945). Hamam, S.E. M., Lu, B. C.-Y., Can. J. Chem. Eng., 52, 238 (1974); J. Chem. Eng. Data, 21, 200 (1976). Harmens. A., Cryogenics, 15, 217 (1975). Horvath, A. L., Chem. fng. Sci., 27, 1185 (1972). Hsi, C., Ph.D. Thesis, University of Ottawa, 1971. Joffe, J., Schroeder. G. M., Zudkevitch, D., AlChE J., 16, 496 (1970). Jones, M. L., Jr., Mage, D. T., Faulkner, R. C., Jr., Katz. D. L., Chem. Eng. Progr. Symp. Ser. No. 44, 59, 61 (1963). Kato, M., Chung, W. K., Lu, B. C.-Y. Chem. Eng. Sci., 31, 733(1976a): Can. J. Chem. Eng., in press (1976b). Lu, B. C.-Y.. Chang, S.-D., Elshayal, I. M., Yu, P., Gravelle, D., Poon, D. P. L., Proceedings, First International Conference on Calorimetry and Thermodynamics, pp 755-766 Warsaw, 1969. Michels, A., Levelt, J. M., De Graff. W., Physica, 24, 659 (1958). Modell, M., Reid, R. C., "Thermodynamics and its Applications", p 535, Prentice-Hail, Englewood Cliffs, N.J., 1974. Parrish, W. R.. Hiza. M. J.. Adv. Cryog. Eng., 19, 300 (1974).
Passut, C. A., Danner, R. P. lnd. Eng. Chem., Process. Des. Dev., 12, 365 (1973). Peng, D.-Y.. Robinson, D. B. lnd. Eng. Chem., Fundam., 15, 59 (1976). Prausnitz, J. M. "Molecular Thermodynamics of Fluid-Phase Equilibria" Prentice-Hall, Englewood Cliffs, N.J., 1969. Prydz. R., Goodwin, R. D., J. Chem. Thermodyn., 4, 126 (1972). Rediich, O., Kwong, J. N. S.,Chem. Rev., 44, 233 (1949). Robinson, D. E., Kaka, H., Proceedings of 53rd Annual Convention of GPA, pp 14-20, Denver, March 1974. Sage, B. H., Lacey, W. N., API Research Project No. 37 "Thermodynamic Properties of the Lighter Paraffin Hydrocarbons and Nitrogen" (1950); "Some Properties of the Lighter Hydrocarbons Hydrogen Sulphide and Carbon Dioxide", (1955). Simonet, R.. Behar. E. Chem. Eng. Sci., 31 37 (1976). Soave. G.. Chem. Eng. Sci., 27, 1197 (1972). Stryjek, R., Chappelear. P. S.,Kobayashi. R., J. Chem. Eng. Data, 19, 334 (1974). Toyama, A.. Chappelear, P. S., Kobayashi, R., Adv. Cryog. Eng., 7, 125 (1962). Vogi, W. F., Hail, K. R.. AlChE J., 16, 1103 (1970). Webber, L. A., NBS J. Res., 74A, 93 (1970). Wichterle, I., Kobayashi, R.. J. Chem. Eng. Data, 17, 4 (1972a); 17, 9 (1972b). Wilson, G. M., Silverberg, P. M., Zeilner. M. G., U.S. Departeent of Commerce, Technical Documentary Report No. APL TDR -64-64, Apr 1964. Wilson, G. M., Adv. Cryog. fng., 11, 392 (1966). Young, S.,Sci. Proc. Royal Dublin SOC., 12, 374, (1910). Zudkevitch, D., Joffe, J., AlChEJ.. 16, 112 (1970).
Received /or reuieub December 8, 1975 Accepted August 4, 1976
A Simple Method for Determining the Reaction Rate Constants of Monomolecular Reaction Systems from Experimental Data Feg-Wen Chang, Thomas J. Fitzgerald," and Jin Yong Park Department of Chemical Engineering, Oregon State University, Corvallis, Oregon 9 733 I
A new method is presented for obtaining reaction rate constants for any set of simultaneous first-order chemical reactions. For k independent reactions the required data consist of measurementson k species at k 1 regularly spaced points in time. The method is exact and gives precise results for error-free data. Redundant data are treated by a simple least-square regression technique.
+
Introduction The monomolecular reaction system is of particular interest to chemical engineers because many important chemical processes are first-order or approximately first-order reactions. Wei and Prater (1962) used the characteristic direction method to determine the reaction rate constants of monomolecular reaction systems. Although the method shown in their work demonstrated a better way for determining the reaction rate constants than the conventional curve-fitting techniques (Haag and Pines, 1960), it is still tedious and complicated, particularly in reaction systems with many chemical species. The object of this paper is to develop a method which provides a more accurate and much easier way for determining the reaction rate constants from experimental data than other previously published methods.
this system is shown in Figure 1.The ith species of this system is designated by A, and the concentration by c,. The reaction rate constant of the ith species to the j t h is denoted by k,,, i.e.
Determination of Reaction Rate Constants from Experimental Data Consider a reversible (or irreversible) monomolecular reaction system of n chemical species. A schematic diagram of
Define
From the mass balance of the j t h species at time t , we find dt
Ind. Eng. Chern., Process Des. Dev., Vol. 16, No. 1, 1977
59
Figure 1. A schematic diagram of a monomolecular reaction system.
Figure 2. A schematic diagram of a monomolecular reaction system for a discrete time basis.
Substituting eq 2 into eq 1,we obtain
(i, j = 1 , 2 , . , . , n )
kl,* = k y d A t
(8)
Substituting eq 8 into eq 7 we get The above equation can be written in matrix notation as
c,(t
+ At) n
where C ( t )is the concentration row vector at time t , i.e., C ( t ) = [ c l ( t )c * ( t ) . . . c , , ( t ) ] ,while K is the reaction rate matrix with elements k,,, i.e.
L Taking the Laplace transform of eq 4 with respect to t , we get
sC(S)- C(0) = C ( s ) K
C(s) (SZ
- K ) = C(0)
C(s) = C(0) (SI - K)-1
+ +
(10)
+
where C ( t A t ) is the concentration row vector at time t At, i.e., C ( t A t ) = [ c l ( t At)c*(t A t ) . . . c,(t A t ) ] ,while Kd is the reaction rate matrix for a discrete time basis with elements k,,*, i.e.
+
+
+
Thus
C ( t -k 2At) = c(t + At)Kd = C ( t ) K d Z
c(t + 3 A t ) = C ( t + 2At)Kd = C(t)Kd3 and in general
C(t
+ mAt)
(for all t 2 0; rn = 0 , 1 , 2 , . . .)
= C(t)Kdm
C ( t )= C(0)eKt
(11)
(for all t 2 0 )
(6)
The reaction rate constant from the ith species to the ith itself, Le., klLd,is nonzero in this case. Let one time interval be At. Then from the mass balance of the j t h species at time t At, we find
+
+ At)
2 (k,,d~t)c,(t)
(for all t 2 0;j = 1 , 2 , . . . , n ) ( 7 )
where c, ( t )is the concentration of the ith species at time t , while c, ( t A t ) is the concentration of the j t h species at time t At. Define
+
Ind. Eng. Chem., Process Des. Dev., Vol. 16, No. 1, 1977
C ( m )=
[ C l h M r n ) .
,
. cn(m)l
Now let us find the reaction rate matrix K of the continuous system that will give the concentrations that are obtained from the matrix Kd for the discrete system at t = 0, A t , 2 A t , . . . . By comparison of eq 6 and 12 when t = mat, we know that pKAt
= Kd
(13)
or
A, kVd_ AI
1=1
(12)
where C ( m ) is concentration row vector a t time t = mAt, i.e.
Now, suppose that we would like to treat the same reaction system on a discrete time basis instead of a continuous time basis. A schematic diagram of this system for a discrete time basis is shown in Figure 2. The reaction rate constant of the ith species to the j t h for a discrete time basis is denoted by k l l d ,i.e.
=
( m = 0 , 1,2 , . . .)
C ( m ) = C)O)Kd"
(5c)
Inverting the Laplace transform of eq 5c, we obtain
60
+
c(t A t ) = C(t)Kd
In a shorthand, eq 11 can also be written as -
+
The above equation can be written in matrix form as
(5b)
and
c,(t
1=1
(54
where c ( s ) is the transformed concentration row vector. Thus
(for all t 2 0; j = 1, 2 , . . . , n ) (9)
k,,*c,(t)
=
(4)
Equation 14 simply tells that if we can find the matrix Kd from which In Kd is computable, then we can find the true matrix K. Next, we develop a method to determine matrix Kd for a discrete time basis. From eq 12, we know (31)
=
('(O)K,,
c'(2) = I'(0)KdL= ('(1)Kd ('(n)=C(O)K,,"= C ( n - 1)Kd
(I&)
FW1= FW&d
Obviously, eq 15a can be written as
Now since FWo is a square matrix, it can be inverted to give Kd =
If the first square matrix in the right-hand side of eq 15b is nonsingular, then the unique solution of eq 15b is
Kd
(19)
Any matrix F can be used so long as FWo is not singular, and if no data errors were present, the precise values for the matrix Kd would be obtained. In practice, however, some F's are better than others. A particularly good (and simple) choice for F is to let F be the transpose of WO.Thus
(16)
=
(FWo)-'(FWl)
Kd = ( WoTWo)-l(WoTW1)
(20)
This choice of F will minimize the mean square of the difference between the measured and predicted concentrations. This can be shown as follows. Let y represent any column in the matrix W1 and let x represent the corresponding column in the matrix Kd. Then y represents the actual concentrations of a species and Wox represents the values predicted for that species. The error is thus ( y - Wox).The mean square error is proportional to
or
Kd =
ZI(Y
- Wox)12= (Y - Wox)T(y- Wox)
Define g ( x ) = (Y - Wox)T(y- Wox)
Hence, eq 17 and 14 simply indicate that in order to find reaction rate constants h,, of a monomolecular reaction system consisting of n species, the only information we need to know is n ( n - 1)transient concentration measurements (possibly including the initial conditions). The initial concentration row vector C(0) must be chosen so that the concentration matrix will not be singular. For instance, choosing zeros as the initial conditions for a pair of simple side-products will produce a singular matrix. Generally, however, the matrix will not be singular, and an easy remedy is to use a different set of initial concentrations. Reducing t h e Effect of Measurement Error Measurement errors occur in the collection of real data. So long as these errors are random (and not systematic), their effect can be reduced by using additional data points. Suppose that the concentrations of all n species are measured at time 0, A t , 2At, . . . M A t (where M > n). The first augmented concentration matrix Wo which is similar to the square matrix in the right-hand side of eq 15b is
w,= C(2)
Taking the derivatives with respect to each element in x gives g ' ( x ) = -2(yTWo
- xTWoTWol
w1= W
8d (18) Premultiplying each side by an arbitrary ( n X M ) matrix F , we obtain
(23)
Setting g ' ( x ) = 0, we have
yTWo - XTWoTWo = 0 Thus X T
= (yTWo)(W07'Wo)-'
or x = (WoTW0)-1(WoTy)
(24)
Since y is any column in W1 and x is the corresponding column in Kd, it follows that the sum of the squares of the errors between the measured and predicted concentrations is minimized when Kd is found according to the formula
Kd = (woTWo)-'(WoTW1)
(20)
A Numerical Example Suppose that a reversible (or irreversible) monomolecular reaction system consisting of three hypothetical chemical species is given by
(25)
Similarly, the second augmented matrix W1 which is similar to the square matrix in the left-hand side of eq 15b is
Thus
(22)
In order to obtain the transient concentrations in this system, we assume the values of reaction rate constants but pretend that only the experimental data obtained from that system are available. Then the reaction rate constants computed from the method presented in this study can be compared to the actual reaction rate constants. In this manner, it is possible to demonstrate how well the method works. For our example, we choose k12 = 0.45; k ~ = 1 0.15; k31 = 0.20; K 19 = 1.20;1223 = 0.60. The principle of detailed balancing (Amundson, 1966) which in effect states that equilibrium is independent of path, requires that we let Ind. Eng. Chem., Process Des. Dev., Vol. 16, No. 1 , 1977
61
Table I. Transient Concentrations of Three Species in the System
The true value of matrix K is
hi3
+ kz)
-(kn
0.5000 0.2560 0.1602 0.1228 0.1084 0.1029 0.1009 0.1002 0.1000 0.1000
0 0.5 1.o
1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
0.5000 0.4328 0.3842 0.3520 0.3316 0.3190 0.3113 0.3067 0.3040 0.3024 0.3014
0.1000 0.1000
-
m
k:iz
0
0.3112 0.4556 0.5252 0.5600 0.5781 0.5877 0.5930 0.5960 0.5977 0.5987
-
k31k12
= k 2 3 -= (0.60) k21k13
0.2oOo
01
C(2) = 10.1084 0.3316
01086
0.3319
0.1619
0.3833
0.4549
. 01086
0.3319
05595
0.1Oll
03111
0.587j
[
0.56001 111
Kd
Substituting into eq 17, we find
I
0.3842
0.1a2
0.3842
0.1084 0.1009
0.3316 0.3113
0.4656 0.5600 0.4556 0.5600
0.3316
I
]
0.2410
0.2260 0.5424 0.1910
0.5281 0.3831 0.7204
We need In Kd to calculate the reaction rate matrix K . From the polynomial method on matrices given in texts (for instance, Bronson, 19691, we know that In Kd can be expressed as a matrix polynomial of degree n - 1
+ alKd + a d +
In X I = a 2 X I 2 alXl
(26)
(27)
where A, are the eigenvalues of matrix Kd. From the characteristic equation of matrix Kd, i.e., det(Kd - XI) = 0, we find XI = 0.9998; Xz = 0.3484; X3 = 0.1586. Substituting the values of X I into eq 27 to determine the a,, we find a2 = -3.0052; a1 = 5.6699; a0 = -2.6650. Substituting the values of a, into eq 26, we obtain 3
0.29%
L203.51 0.5970 -0.4991
Substituting this into eq 14, we get
[
-1.6515
K=
62
01516 0.1992
-1.56% 01016 0.2117
0.4442 -0.7456 0-2985
0.5410
05656
0.3680
0.1815
0.7271
03495 -0.6763 0.2762
0.3495
-0.6763 0.2762
L2198 0.5695 -0.4869
z:yi
W" = ~
C( C (5)6 4
0.1592
LZO
0.5272
0.3064
05560
0.4 318
03122
0.1592 01206 01078 0.1049 01011
0.3860
0.4548
0322
0.5212
0.3322
0.5600 05767 0.5873 05931 0.5951
0.1049
C(7)
-
03184 0.3116
05970
0.1w
-0.4991
_0.0988
0.3055
Ind. Eng. Chem., Process Des. Dev., Vol. 16, No. 1, 1977
0.4548
0.5600 05767 0.5873 0.5931-
i
and the matrix W1 is
0.3860
0.3184 0.3116 0.3064
120351
1 1
L2198 0.5695 -0.4869
-
C(2)
+ a0
0.4442 -0.7456
1
om10
By comparison of (K)calcand (K)true,we find that the agreement is reasonably good. The largest relative error is about 30% (for k Z 1 ,which is the smallest rate constant). Case 3. Data with E r r o r , Least-Square Fit for Redundant Measurements. Here, we use five extra sets of data and let A t = 0.5. 'ILen the matrix W Ois r 3 O m 3 0.0004 0.2560 0.4318 03122
where the a, 's are given by the equation
-1.6515 0.l516 0.1992
[
K=
00764 [0:@77
0.5877
In Kd = azKd2
[l 0.0633
Substituting this into eq 14, we get
I
0.1M2 0.1084
=
-1.5635 01016 0.2117
=
C(3) = [0.1009 0.3113 0.58771
K,,=
0.5595
Then using the same computation procedures as Case 1,we obtain
C(1) = [0.1602 0.3842 0.45561
0
L2000] 0.6000 -0.5000
(0.20)(0.45) = 0.3 (0.15)(1.20)
C(0) = [0.5000 0.5000
0.5000
0.4500
-0.7500 0.3000
The agreement of course is good; the differences occur only because of limited precision in the calculations. Case 2. Data with Error, No Redundant Measurements. We now add a random error with a standard deviation of 0.25% maximum amplitude and correct the sum of each row to one. Again, let At = 1.0. Substituting the noisy data into eq 17, we find
A number of transient concentrations of these three species assumed to be measured from the experiment are tabulated in Table I. Case 1. T h e D a t a A r e Error-Free. Let At = 1.0. Then from Table I, we know
O . m
+k d -1.6500 0-1500
0.6000
0.3000
kz
k32
-
-
C(s) = row vector of the Laplace transform of concentra-
Substituting them into eq 20, we find 0.4710
01500
O.W]
K,j = ( W O ~ W , ) - ' ( W , , ~ W , )0.0410
0.7200
0.2400
0.1160
0.8160
[O.O6W
Then using the same computation procedures as Case 1,we obtain 0.5893
1
-0.7976
OX35
0.0531
-0.3551
0.3051
0.1075
0.1451
-0.2489
Substituting this into eq 14, we get -1.5952
K = - _Ind K= - ! In K ilt 0.5
=
01062 0.2133
0.42'70
1.1786
-0.7102 0.2902
-0.4958
0.61.12
1
As expected, we find that the agreement is better by using extra data points. The relative error for the smallest coefficient is still about 30% but the other coefficients agree quite well. Nomenclature
A , = the ith species in the reaction system a, = constant in eq 26 and 27 c, ( t ) = concentration of the i t h species a t time t c , ( m ) = concentration of the ith species a t time mAt C ( t ) = concentration row vector a t time t C ( m ) = concentration row vector a t time m A t
tion
F = an arbitrary matrix g = function I = identitymatrix ki; = reaction rate constant of the i t h species to the j t h kijd = reaction rate constant of the ith species to the j t h for a discrete time basis hi;* = defined byeq 8 K = reaction rate matrix Kd = reaction rate matrix for a discrete time basis n = number of chemical species s = Laplacevariable t = time W O = augmented concentration matrix W1 = augmented concentration matrix x = column vector in the matrix Kd y = column vector in the matrix W1
Greek Letters At = one time interval X = eigenvalue of matrix Kd Literature Cited Amundson, N. R., "Mathematical Methods in Chemical Engineering-Matrices and Their Application," p 230, Prentice-Hall, Engiewood Cliffs, N.Y., 1966. Bronson, R., "Matrix Methods-An Introduction," p 110, Academic Press, New York, N.Y., 1969. Haag, W. O., Pines, H., J. Am. Chem. Soc., 82, 387, 2488 (1960). Wei, J., Prater, C. D., Adv. Catal., 13, 203 (1962).
Received tor revieu December 4, 1975 Accepted .July 23, 1976
Design Correlations for Autothermal Reactors with Internal Countercurrent Heat Exchange Jaime P. Ampaya and Robert G. Rinker' Chemical and Nuclear Engineering Department, University of California, Santa Barbara, Santa Barbara, California 93 106
Correlations are developed to relate the critical feed temperature, corresponding to blow-off conditions, and the critical catalyst-bed entrance temperature to the operating and design parameters of autothermal reactors with internal countercurrent heat exchange. The proximity of the blow-off conditions to the optimum point of operation makes the correlation of the critical operating temperatures highly significant. Use is made of both irreversible and reversible first-order reactions. The capabilities of the correlations are illustrated with the water-gas shift reaction treated as pseudo first order.
Introduction In designing chemical reactors, removal of the heat of reaction is often desirable to limit the reaction temperature, to recover the energy produced, and for an equilibrium reaction, to minimize the effects of the reverse reaction. The last reason is probably the most important because of its relationship to the optimum temperature profile needed for maximum production rate. Attempts a t simulating the optimum profile, which may be too costly to actually attain in practice, have used inter-bed cooling as in converters for synthesis of sulfuric acid and have used injection of cold feed as in water-gas shift reactors. Another approach and indeed a useful one is to use an autothermal reactor with internal countercurrent heat exchange, referred to hereafter as an ARICHE. Perhaps the
most important industrial reactor of this type is the HaberBosch ammonia converter. The feed, serving as the coolant to an ARICHE, runs countercurrent to the effluent stream and is preheated by the exothermic reaction. The results of calculations based on a one-dimensional, plug-flow model of an ARICHE as well as a conventional adiabatic reactor under conditions of identical reactor size, space velocity, and reaction parameters are shown in Figure 1. I t can be observed that (1)the maximum attainable conversion for an ARICHE is higher than that for an adiabatic reactor; (2) for any particular feed composition, there exists an optimum reactor feed temperature for producing a maximum conversion in an ARICHE; (3) conversion for an ARICHE can be multi-valued; and (4) an ARICHE has a minimum or critical Ind. Eng. Chem., Process Des. Dev., Vol. 16, No. 1, 1977
63
