ASYMPTOTIC INTEGRATION OF OVER-ALL MASS AND ENERGY BALANCES FOR
NONISOTHERMAL REACTION SYSTEMS O W E N T. H A N N A ’ A N D Rensselaer Polytechnic Institute, Troy, N .
ROBERT S.
K A P N E R
Y.
A technique is presented for integrating over-all mass and energy balances associated with certain nonisothermal reaction systems. In the case of adiabatic systems, the technique yields results which are asymptotically valid for small times and also for large values of activation energy at all times; energies of value no greater than those found typically permit good results. The method can be used for complicated kinetics and variable thermal properties and as such constitutes the only analytical method available for most adiabatic reactions. It also can be used for nonadiabatic systems, if the heat transfer rate is sufficiently small.
HE nonlinear Arrhenius expression, k = k , exp(-E/RT), Tis a major factor prohibiting analytical solution of many nonisothermal reactor design problems, particularly for adiabatic reactor design as well as for reaction systems involving heat transfer across system boundaries but still imposing nonisothermal conditions. Such systems usually require the simultaneous solution of mass and energy balances by numerical means. I t is desirable to obtain simple analytical solutions to such problems in order to determine the dependence of the solution on various reaction parameters. This would permit rapid evaluation of alternative reaction schemes. This paper develops a new method of arriving at approximate solutions to nonisothermal, irreversible reactor design problems of considerable complexity. Recently some attention has been paid to this problem for certain simple cases. One scheme, proposed by Barkelew, relies on approximating the Arrhenius equation by k = c e b ( T - T o ) . This approximation permits a simple integration of some design equations (2, 4 ) . For adiabatic reactions, Aris has shown an exact integration in terms of tabulated exponential integrals ( I ) . Both methods, however, suffer from restrictions which limit their usefulness. The Barkelew approximation can be used for moderate temperature changes, but becomes inaccurate when the temperature range is very large. The exponential integral solution for adiabatic reactions requires that the design equation integrated be of rather simple form. I n cases of large temperature changes and complicated design functions, particularly those involving temperature-dependent thermal properties, these methods cannot be used. The new method described here can, in general, be used for these difficult cases. I n principle, the method depends upon techniques for integrating functions containing large parameters. I n the cases considered the large parameter is E/R, which arises directly from the Arrhenius equation.
Batch Reactor
To focus attention on the solution method we treat first the simplest possible nontrivial system, that of an irreversible, 1
Present address, Airplane Division, Boeing Co., Renton,
Wash. 116
I&EC FUNDAMENTALS
-
constant density batch reaction. For this case, the reaction B products has the rate equation
where the reaction velocity constant is given by the Arrhenius expression k = k, exp(--E/RT)
(2)
A differential energy balance for this system, assuming adiabatic operation and a constant heat of reaction and heat capacity, is
(3) The relationship between reaction time and conversion, xB, is given by the integral of Equation 1
“JJ0
dX B
(4)
The method proposed here for evaluating Equation 4 makes use of the idea that the major contribution to the integral comes from the vicinity of the minimum temperature, since the exponential will be maximum when the temperature is minimum. The procedure is to introduce a new variable of integration which is zero near the lowest temperature and which causes eaIT to become e-au. For example, when the reaction is exothermic, let 1 / T = l/To - u. Then, Equation 4 becomes
. ,
The integral in Equation 5 is of a form which can be represented asymptotically for large values of a by means of Watson’s lemma (6). For our purposes the lemma must be modified slightly for best results. The idea of the lemma is that when a is large, the negative exponential reduces the integrand so much at larger values of u that the value of the integral from zero to any fixed upper limit can be obtained accurately by expanding the non-
exponential part of the integrand in a Maclaurin series and integrating term by term from zero to infinity. This process generates an expansion which usually diverges if interpreted as an infinite series, but has the property that for any fixed number of terms the expansion approximates the original integral to any desired degree of accuracy for sufficiently large values of a. Since the value of a is fixed (and hence cannot be increased so as to achieve better accuracy), the procedure is limited in its usefulness to cases where a happens to be large enough so that the method gives good results. [In situations where the value of a is only moderate, acceptable accuracy may sometimes be achieved by this method when Euler’s transformation is employed ( 3 ) . ] In many kinetics problems where the temperature dependence of the rate constant is important, a simple one- or two-term application of the above procedure, slightly modified, yields very accurate results. Also, it is easy to obtain expressions for the upper bound on the error made by using this technique. Accuracy is normally much better than this upper bound indicates. The procedure yields results which become exact both for small reaction times for any value of a and for large values of a a t any time. T o illustrate the procedure, we first write Equation 5 as
- l/T,),
where for the exothermic case A = l / T o ,R = ( l / T o and
This notation permits us to obtain results valid for both exothermic and endothermic systems, as will be seen. Now expand G(u) in a two-term Taylor expansion and integrate to obtain
t =
[F
(1
- e-aR) + G’(.$)
SOR
~e-‘~du]
(7)
where 5 is some point on the interval of integration. Thus the absolute difference between the exact value o f t and the first approximation to it will be less than the maximum value of the last term on the right of Equation 7 . It is therefore easily shown that
Inequality 8 holds regardless of whether the last factor on the
right is R2 or l / a 2 . Since G’(.$)m.x (shown in Table I) is independent of R and a, Inequality 8 shows that the first approximation will be asymptotically correct for either small R or large a. Since small R implies small ( T , To),the approximation is seen to hold at any a at small times or for any times at large a. In general, further terms of the expansion will improve the accuracy up to a point; however, for simplicity, and because the numerical results do not seem to require it for the examples given, only one or two terms are presented here. T h e three-term result [obtained by replacing G’(E) by G’(0) in Equation 7 and integrating] is given by
-
where
The asymptotic character of the three-term expression can be established as for the one-term result. Equation 9 applies to both endothermic and exothermic reactions; appropriate values of the quantities A, R, G(O), G’(O), and G”(0) are shown in Table I. The integral in Equation 6 can be evaluated for large a according to Watson’s lemma by expanding G(u) and integrating term by term from zero to infinity. This procedure is useful for a fixed value of R, but since here we are concerned also with small values of R (at short times), we do not neglect the e-’R terms in Equation 7 . Thus if the eFaR terms in Equation 7 were set equal to zero, the result would be the first two terms of an asymptotic expansion for large a. However, neglecting the e-aR terms would result in error for R * 0 ; the retention of the e-aR terms eliminates this difficulty. For the example of an adiabatic exothermic reaction with constant C, AH, and density a particular set of numerical values, listed in Table 11, was selected. Using these values, which were chosen to be reasonably typical of many such systems, a numerical integration of the mass and energy balances was performed in order to assess the accuracy of the proposed method. Results of this integration are compared in Table I1 with the predictions of Equation 9 for one, two, and three terms, respectively. Agreement is seen to be good. By reversing the sign of AH while retaining the other values of the parameters, a similar integration and comparison with Equation 9 was made for endothermic conditions, For the
Terms in Asymptotic Solution of Endothermic and Exothermic Reaction Systems
Exothermic l/To
min u = 0 corresponds to final conditions.
u = 0 corresponds to initial conditions,
VOL. 6
NO, 1
FEBRUARY 1967
117
endothermic case the agreement was slightly better than for the exothermic results shown in Table 11.
perature-dependent heat capacities of system components. The design equation for a tubular reactor assuming ideal plug flow is
F l o w Reactor, a Variable-Property Example
Vr=FJ
As a second problem we treat a flow reactor design of considerable complexity involving variable density and temperature-dependent thermal properties. An example of this type of problem, which involves the first-order decomposition of phosphine, is developed in detail by Hougen and Watson (5) and shown in Table 111. Consider the gas phase, irreversible, nth order, decomposition of species B into products A and C, PBk)
-t
OAk)
xB
dXB T,
y(xB,
Substituting 11 and 1 2 into 14, and assuming constant pressure operation
t-1
( 1 0)
f Yc(8)
for which the variable density rate equation, assuming perfect gas behavior, is
dXB
(15)
The usual procedure is to integrate Equation 15 numerically using pairs of values of xB and T from the energy balance, Equation 13. Using the method described in this paper we first transform Equation 15 into
+ -
Here xB is the conversion level of B and 6 = (a y p)/p is the gas phase expansion factor. The reaction velocity constant, which includes a temperature-dependent preexponential factor, is of the form k = (DTb)e-aiT
(12)
where, for this endothermic reaction, A = l / T f ,R = l / T f
- l / T o ,u
= 1/T,
- 1 / T , and
For adiabatic operation a unique relationship between xB and Tis found from an energy balance and may be written XB
=f(T)
( 1 3)
where f ( T ) is a ratio of polynomials in T arising from tem-
The value of the integral in Equation 16 is given by Equation 9. Using only the first term of the expansion of G(u), the oneterm asymptotic representation of the design Equation 15 is
F
v' = DAdiabatic Batch Reactor, Exothermic Reaction Reaction Time, t, Sec. T, Eq. 9 Eq. 9 Eq. 9 ' K. Numerical ( 7 term) (2 terms) ( 3 terms)
Table II.
xB 0.01
602 610 620 660 700 734
0.06548 0.06494 0.06547 0.06548 0.2696 0.2595 0.2692 0.2696 0.10 0.433 0.404 0.431 0.433 0.30 0.667 0.578 0.650 0.663 0.50 0.713 0.601 0.685 0.704 0.67 0.724 0.604 0.691 0.712 Reaction. B + products k, = 1014sec.-l a = E / R = 20,500' K . -AH/C, = 200' K. To = 600" K . ~ ( x B = ) 1 XB 0.05
-
Table 111.
XB
- e-aR)
"eaAG(0)(l a
(1 7)
where
Table I11 lists the numerical values of reaction parameters for this problem. Also shown is the comparison of reactor volumes as found by numerical integration of Equation 15 and as evaluated with the asymptotic result, Equation 17. The agreement is very good for the conversion range investiga ted. Nonadiabatic Problems
Adiabatic F l o w Reactor, Endothermic Reaction Reactor Volume, V,, Cu. Ft. Eq. 9 Eq. 9 T , ' K. Numerical ( 7 term) (2 terms)
7.031 6.988 6.988 0.0117 938.02 70.85 69.30 69.30 0.0531 913.24 40.09(101) 41.32(101) 40.07(101) 0.1005 884.92 27.94(102) 28.88(102) 27.92(102) 0.1539 853.06 17 .29(103) 17 .87(103) 17.27(103) 0.2015 824.74 12.82(105) 13.24(106) 12.81(106) 0.3034 764.56 Reaction (5). 4 PH,,,) + Pqg) 6 Hqg) F = 14.65 lb. moles PHa/hr. P = 2 atm. D = 4.855(1O16)(hr.-' K.-b) a = 43,664' K. 6 = 2 n B o = 1 mole B/mole fed Reaction order = n = 1 6 = 314 T o = 945.1' K. Energy balance. X B = 9145 - 6.70 T - 0.00315 T 2 5.095 T - 0.0019 TS 0.14 (10-e)T3f 4341
+
O
+
118
R (-) P
I&EC FUNDAMENTALS
The foregoing results apply only for adiabatic operation, I t is also possible to consider certain more general nonisothermal problems using the preceding technique. I n the case of batch operation, the energy balance, Equation 3, becomes
- -_
dT dt
T ) I QA' CP NCP
mr(xB,
(18)
Combining Equations 1, 2, and 18 yields the following expressions, which may be used to calculate reaction times, provided that dxB/dT is well behaved.
In nonadiabatic situations, the results given earlier do not generally apply, because the integrand in Equation 19 does not satisfy appropriate conditions. In particular, the quantity l/k$(x,) - dxB/dT is of exponential order. Further, in cases where the temperature goes through a maximum or minimum, the quantity dxB/dT is unbounded. For cases where dT/dxB does not pass through zero, several different asymptotic results can be obtained. Of these results, the greatest practical interest corresponds to problems where the main temperature change is caused by the chemical reaction as opposed to heat transfer. In this connection the problems of exothermic reaction with cooling and endothermic reaction with heating are most important. When heat transfer is of the same order of magnitude as chemical reaction effects-in the case of a temperature maximum, for example-the results given do not apply. The case where heat transfer controls the temperature can be treated, but seems unimportant practically. The following development is based on the assumption that the chemical reaction is the primary source of temperature change and that external heat transfer is of secondary importance. For this case, where heat transfer is taken to represent a perturbation to the adiabatic problem, immediate extension of the adiabatic results can be made. If AH is sufficiently large in magnitude compared to the variable u may be introduced for T again and the expansion of G(u) = I / [ k $ ( x ~ ) ] X (dxB/dT)carried out together with term by term integration. This result, with the exception that G(O), G'(O), and G"(0) now have different values, is the same as Equation 9. The solution so obtained represents a perturbation of the adiabatic problem which is useful for cases where the temperature change is mainly governed by the reaction, and where either a is sufficiently large or ( l / T o - l/T,) is sufficiently small. One difficulty concerning the direct use of Equation 9 is that the final temperature is no longer a simple function of the final composition. For exothermic problems, where initial conditions mainly determine reaction times or volumes, this difficulty is not critical; except for small times or volumes the effect of T , is only secondary. However, even for the exothermic case neglect of the effect of T , results in the loss of some accuracy. O n the oth.er hand, the role of T , in endothermic problems is vital, since here quantities are evaluated a t final conditions. T o determine T , it is generally necessary to solve the conservation equations numerically. We note, however, that Equations 1 , 2, and 18 may be combined (for constant properties) to give
quantity To - AH/Cp(xl - xo), which is obtained by neglecting the contribution of the integral for small R and 4. This procedure yields accurate values of T f for small conversions. For larger conversions, the initial approximation to Ti substituted into the integral in Equation 22 will be inaccurate; however, in this case the value of R is larger, and the asymptotic representation of the integral is sensitive to R only when R is small. The applicability of these results is shown in Figures 1 and 2. T h e results shown there correspond to the constant property, adiabatic system referred to in Table 11, except that now 0 = h(T, - T ) . The figures indicate that for this I
'
I
I
I
I
I
ASYMPTOTIC - - - - e - - - NUMERICAL
/i R
Xf= 0.5
/
4,
I
e
0.9
I
I
L
2
0.5
W
a
0.1
'
0
___--*
---
,-
-,Xt=O.l
I
I
I
I
I
04
00
12
16
20
I
HEAT TRANSFER COEFFICIENT ( h )
Figure 1 . Numerical and asymptotic calculation of batch reaction time in exothermic reaction with cooling
I
I
I
I
ASYMPTOTIC - - - e NUMERICAL
\
6801 W LT
3
I
I
I
--
----. Xf=0.5
n I-
Using the previous notation, we have
I W
I-
z
0 t
0
OG(u)e-audu
(22)
a W
a
Xf'O.l
- 026
T h e integral in Equation 22 is of exactly the same form as we have treated previously, provided we take the product QG(u) in Equation 22 to be G(u) in Equation 9. Thus T , can be obtained from the asymptotic expansion of the integral in Equation 22. The quantity R depends on T,; it is convenient to replace the value of T , occurring in R by the
0
0.4
0.8
1.2
1.6
20
HEAT TRANSFER COEFFICIENT ( h )
Figure 2. Numerical and asymptotic calculation of batch temperature in exothermic reaction with cooling VOL. 6
NO. 1 F E B R U A R Y 1 9 6 7 119
example useful results are achieved for h as large as 2 . Apparently the method used here will be less accurate for endothermic reactions with heating, since here the residence time is sensitive to the value of T,. T o reiterate, the preceding technique applies for sufficiently large and small values of a and respectively, as well as for any a and 4 at sufficiently small temperature changes (conversions). Thus, generally speaking, for larger conversions the above method is best suited to problems involving steep temperature variations. For example, many problems involving exothermic reaction with heating or endothermic reaction with cooling, although of less practical interest, could be handled nicely with the above procedure. The example shown in Figures 1 and 2 , although less suited to the requirements of the technique, represents a problem of greater practical interest. With regard to estimating the error incurred by using a certain number of terms in Equation 9, a general heuristic rule is that the error is of the same order of magnitude as the first neglected term-that is, if the first neglected term is not somewhat less than the sum of the preceding terms, the results should be vieLved with suspicion. I n these cases it is sometimes useful to employ Euler’s transformation ( 3 ) .
a,
Conclusions
The integration method proposed here for the purpose of handling nonisothermal reaction problems offers a number of important advantages. First, it is not restricted to simple kinetics or to constant C, AH, or density, which makes it the only analytical method available for most nonisothermal problems. As shown in the examples above, the accuracy obtainable for typical problems is very good. Also, in many cases, determination of the upper bound of the error is not difficult. The method should be very useful where quick estimates are required, especially for complex reaction schemes. The results are particularly convenient for computer work, since the proposed relations will lend themselves to parametric studies and to dynamic control where a simple and accurate characterization of a reaction process is important. Even when detailed numerical integrations of the over-all mass and energy balances are performed, the present results should be valuable for checking and interpolation. Also, advantage
could be taken of this technique for experimental kinetic studies involving integral reactors. Nomenclature a
= E/R
A’ A
= =
heat transfer area 1/T,or l / T f b = constant in Arrhenius Equation 12 C, = heat capacity 6 = expansion factor D = pre-exponential factor in Arrhenius Equation 12 E = Arrhenius equation activation energy F = feed rate h = convective heat transfer coefficient AH = heat of reaction (-for exothermic reaction) k = reaction velocity constant n = reaction order nB = moles of component B/mass fed N = moles charged to batch reactor P = pressure 4 = convective heat flux to batch reactor r = reaction rate equation R = gas constant, or upper limit of integral t = reaction time T = absolute temperature T, = wall temperature V, = reactor volume XB = conversion, moles B converted/moles fed u = l / T o - 1 / T (exothermic), l / T f - 1 / T (endothermic)
SUBSCRIPTS o
f
= =
initial condition final condition
literature Cited
(1) Aris, R., “Introduction to the Analysis of Chemical Reactors,” p. 234, Prentice-Hall, Englewood Cliffs, N. J., 1965. ( 2 ) Barkelew, C. H., Chem. Eng. Progr. Symp. Ser. 5 5 , No. 25, 37 (1 959). (3) Bromwich, T. J., Macrobert, T. M., “Introduction to the Theory of Infinite Series,” Macmillan, London, 1949. (4) Groves, F. R., Jr., IND. ENG. CHEM.FUNDAMENTALS 4, 98 (1965). (5) Hougen, 0. A., Watson, K. M., “Chemical Process Principles,” 1st ed., Vol. 111, pp. 865-9, Wiley, New York, 1947. (6) Jeffreys, H. J., Jeffreys, B. S., “Methods of Mathematical Physics,” p. 501, Cambridge, New York, 1946. RECEIVED for review May 28, 1965 ACCEPTEDSeptember 20, 1966
TRANSIENT RESPONSE AND MOMENTS ANALYSIS OF BACKFLOW CELL MODEL FOR FLOW SYSTEMS WITH LONGITUDINAL MIXING M . H . R O E M E R I A N D L.D. D U R B I N
Department of Chemical Engineering, Texas A&M University, Collexe Station, Tex.
or axial mixing effects in agitated flow procare generally characterized by a differentially continuous diffusion model with the degree of mixing determined by the dispersion coefficient or “eddy diffusivity,” D. ONOITUDINAL
L esses 1
120
Present address, Celanese Corp. of America, Clarkwood, Tex. I&EC FUNDAMENTALS
The relative degree of mixing between two systems may be characterized by the dimensionless Peclet ratio, Pe = (UL/D), of (velocity X length) to eddy diffusivity. A considerable amount of work concerned with theoretical and experimental analyses based on the diffusion model has been published within the past decade (7,3,4, 7,8, 70, 72, 78, 79). The simple serial s tined cell model (9, 74) ha,s also been applied to mixing
