ENGINEERING, DESIGN, AND PROCESS DEVELOPMENT Figure 2 shows the agreement of this relationship with the experimental data. Although in this investigation the roll dimensions were not varied, the complete agreement in form between this correlation and that derived from two roll studies indicates that the conclusions drawn from either of these studies might be applied to the other. I n the two roll study, roll dimensions were varied, and the variation can be expressed according to Equations 1 and 2. I n this study, the density was varied and the variation can be expressed according to Equation 6. Since the equations are the same in form and are a result of derivation and experimental confirmation, conclusions drawn from one case can be applied t o the second case. The difference in the constant of proportionality seems to be directly related to mill loading alone and in no way serves to describe the hydrodynamics of the roll mill beyond mill loading. Because of the large number of quantities measured in this work, the cumulative experimental error was quite large. Among the most important of these errors was the measurement of the total clearance where the problems were the same as those described in the first paper ( 8 ) . Another point a t which an important error may have developed was in the determination of the power consumed by the material on the rolls from the power input to the motor. There was also a question as to whether the temperatures measured to get the viscosity of the material are actually close to those prevailing in the nips where the viscosity is of greatest importance. The correlations obtained are considered to be very good in view of these possible errors. Conclusions
A 4 X 8 inch, three roll, laboratory printing ink mill with floating center roll was studied in an effort to correlate the operating variables, power consumption, mill speed, roll speed ratio, nip clearance, and viscosity. Roll dimensions and density were included in the correlation but not varied in this investigation.
A general correlation equation was derived by summation of separate equations for the feed and apron nips alone which were developed in a previous paper. This equation is
+ 1)1'S(R + k)1'3
where K = ( R ) ( R
Best fit with the measured data was obtained when
+ 1)1/3(R+ k)Ij3
K = (R)(R
This study, in confirmation of earlier. work, has shown t h a t three roll mill operation may be treated as a fluid flow phenomenon. The principles of fluid dynamics that have been employed with great success in the study of numerous chemical engineering operations may now be applied to the improvement of roll mill design and operation. Acknowledgment
The authors wish to acknowledge the combined support of the National Printing Ink Research Institute and the Troy Engine and Machine Co. which made this work possible. literature cited (1) Lower, G. W., Walker, W. C., and Zettlemoyer, A. C., J. Colloid Sci., 8 (l),116-29 (1953). (2) Maus, L., Walker, W. C., and Zettlemoyer, A. C., IND.ENG. CHEM.,47, 696 (1955). (3) Maus, L., Zettlemoyer, A. C., and Gamble, E., I b i d . , 701 (1955). RECEIVED for review December 16, 1953. ACCEPTED October 6, 1954. From a dissertation submitted by James H. Taylor to the Graduate School of Lehigh University in partial fulfillment of the requirements for the degree of Master of Science, June 1952.
Design Equations f o r .
Continuous Stirred-Tan k Reactors Transients of Second-Order Chemical Reactions F. s. ACTON
AND
L. LAPIDUS
Princeton University, Princeton,
I
N. J.
N RECENT years the development of design equations describing the behavior of a chemical reaction in a series of continuous homogeneous stirred-tank reactors has reached a status of maturity. Since the contributions of Denbigh ( 1 , 2 ) both the steady-state operation and the transient approach to the steady state have been considered ( 4 , 6-8, 1 2 ) . Because of the algebraic nature of the steady-state equations it has been shown by Eldridge and Piret (3') and Jones ( 5 ) that the effluent from any reactor may be calculated independently of the order of the chemical reaction. By combining graphical and numeriEal techniques, the system of equations can always be obtained. While experimental data to check the theoretical equations have been few, the agreement between theory and experiment has been excellent. The transient approach to steady state, a condition of utmost physical importance, can be treated by present methods only when certain simplified chemical reactions are occurring (9). For any reaction other than first order or those equivalent to first order, the differential equations become nonlinear. The
706
case of a second-order chemical reaction has only been touched and no solution has appeared in literature. This problem is not unique to stirred reactors. I n many chemical engineering operations the transition from a first-order reaction to one of higher order results in nontractable differential equations (11). The equations presented in this article aid in the design of systems involving transient continuous stirred-tank reactors with a second-order chemical reaction. Solutions of the proposed equations are presented in closed form for the first reactor. For the reactors beyond the first, a set of approximate equations are proposed that allow the effluent concentrations to be calculated within the limits of engineering accuracy. Contrary to the current trend of using an electronic computer to evaluate all interesting combinations of the physical parameters, the computer is used in this work as an intermediate tool t o suggest analytical approximations. This technique deserves emphasis because of the saving in time and money and the more compact form of the results.
INDUSTRIAL AND ENGINEERING CHEMISTRY
Vol. 47, No. 4
ENGINEERING, DESIGN, AND PROCESS DEVELOPMENT Design equations are developed for prediction of transient effluent concentration from any reactor
Consider a series of continuous stirred reactors in which a homogeneous liquid phase second-order chemidal reaction of the following form is occurring
A
+ B (in excess)
--t
products
T h e feed to the first reactor is assumed constant and contains reactants A and B a t concentrations a. and bo, respectively. If complete and instantaneous mixing of the influent t o each reactor is assumed, the effluent from the nth vessel, as well as the feed 1 vessel, will be equal t o the concentration of material to the n in the nth vessel. A material balance for reactant A around the nth reactor during the approach t o steady-state operation results in the equation
+
h consistent set of units must be employed. One may remove b, from Equation 1 by defining E as the absolute excess of B above A . E will be constant from any reactor, and thus E = bo - ao
bn
-
(2)
an
Equation 1 becomes pan-] - pa, =
+
1' 1 da 'kVa,,(E
+ an)
dt
be solved by interpolation. This philosophy unfortunately founders on the stubborn rock of three parameters. T o explore adequately the physical ranges of two parameters is usually feasible, but three can easily require large quantities of expensive computations and lead to a whole book of graphs which no editor would publish. Nor would engineers really care to have this book cluttering up their desks. Accordingly, numerical integration was used to examine a few cases-far too few to allow interpolation between them-in the hope that they would suggest adequate approximations. I n many engineering problems, numerical methods are employed t o solve specific and otherwise apparently untractable equations. Thus far, however, very little exploratory use seems to have been made of computation in order to find analytical approximations for these solutionsapproximations t h a t take up less space than a catalog and are more useful than the few isolated solutions which suggested them. The rather considerable leeway allowed by the tentative nature of many engineering problems makes numerical searches for analytical engineering solutions both attractive and profitable. The present example is offered in support of this thesis. Summary of Results. The equations for the effluent concentration from each of a series of vessels in which a second-order reaction is taking place have been given as Equations 3 and 4. where the vessels are assumed to contain solvents initially but no reactants, and the input concentration of the reactants to the first vessel is kept constant. For computational expediency it i.s desirable to regroup the parameters by defining
or, on rearranging and setting V / q equal to 0
K = ke B = l + K E
T
t/e
=
so that Equations 3 and 4 become
I n addition, the initial condition may be specified by
dadt ,(T)+ ~ a : +, ~a~
a,(0) = b,(O) = 0
(4)
= an-l(~)
n = I, 2,
(6)
an = constant Equations 3 and 4 represent the system of equations that is considered in this paper. Equation 3 is a nonlinear differentialdifference equation that has defined solution in closed form for t h e nth reactor. For the first reactor, however, a closed solution is possible, with Equation 3 becoming (3')
where
2a0 coth (Qot/ae)
+ keE + Qo
Q~ = [(I
+ k0El2 + 4keaoi1/2
U"(0) = 0
a0
=
constant
I n this notation, the exact solution for the first vessel (n = 1 ) is
(7) where The approximate solutions t o Equation 6 use the equilibrium or steady-state (2' = m ) concentrations in their formulas. These are found recursively as
The solution to this Ricatti equation is = 1
. .
15)
Equation 5 may be obtained by standard techniques. The term al(t/e) describes the effluent concentration of reactant A from reactor 1. When t becomes large, this solution approaches the steady-state equation in the literature. The same approach to the second reactor equation is not useful since it contains ul(t)as a coefficient. al(t)is a variable function of time and is sufficiently complicated to preclude closed form solutions. The same sort of remark applies to all subsequent reactors as well, and hence approximate numerical techniques are needed.
&.E
=
where
2an-1,E B Qn-1
+
Qn-1
=
n = 1, 2,
(B'
4KUm-i.~)"'
and are the successive solutions of the quadratic equations obtained by setting
nnequal to zero in Equation 6. dU
By using the definitions of this section, the approximate solutions to Equation 6 for small values of T are
Numerical integration i s employed to suggest adequate approximation
Since the solution of Equation 3 is not possible except for the first reactor, one obvious expedient is t o integrate several typica! sets of the physical parameters, E, Ice, and UO, numerically. If t h e practically interesting values of the pa.rameter space could be spanned, a catalog of curves would then enable any problem to April 1955
where S = BT.
I ND U ST R I A L A ND EN G IN EE R IN G CHE M I STR Y
707
ENGINEERING, DESIGN, AND PROCESS DEVELOPMENT
da, + Ka;: + Ba,
dT
= an-l
'n = 1, 2,
.. .
(11)
where K = 8k B = 1 f K E T = t / O = time measured in units of the tank holdup time When n is one, this equation can be solved explicitly by rewriting it as
where z is a dummy variable of integration. When uois a known constant, this integral is in all standard tables of integrals.. Substituting the two limits, integrating and simplifying,
Figure 1. 1.
1
2.
1
3.
1
Variation of functions for different values of x
- e-z - ( 1 + x ) e-" - - 2 ~ e -~ e-2"
5.
1
- (2
+ x2)e-2 + e-2z
When this equation is solved for ul, Equation 7 is obtained. Since the solution of the problem is known for the steady state, and since the initial conditions are all zero, the problem was t o find an interpolatory function which grows approximately as does the solution t o Equation 3 and has the same limiting values. Two general forms were suggested by limiting cases of Equation 3-the first arising whenever the value of a, is small enough so that the quadratic term, Ka;, is considerably smaller than Ban; the second arising when no reaction takes place (K=O) and only mixing dilutions occur. The equilibrium (steady-state) concentrations are found by . da, setting d~ equal to zero and solving the resulting system of quadratics.
2am-1,E = ___B Qn-I
+
U7i.E
a0.E
where
For larger values of T,the approximate solutions of Equation 6 are a l ( T ) = a l ~ [l
e-Ro]
a,(?") = a 2 ~ [ 1 - (1
+ Rl)e-R~]
(10)
where R, = QnT.
Derivation of approximations
For mathematical manipulation it is more convenient to rewrite Equation 3 as
708
(B2
+ 4Kan-1.~)"2
(The usual formula for the solution of a quadratic gives a computationally annoying form involving the difference between two similar quantities, B and Q . The form given by Equation 12 uses only the sum and is therefore satisfactory for slide rule computation.) When K is zero, B is one, and the sequence of equations daw @
The rational for deriving these approximations is presented in the following section. Both Equations 9 and 10 use identical exponential functions of a single parameter. These are displayed graphically in Figure 1. The computations for the cycle history of a vessel system would be begun by finding the steady-state values from Equation 8 and saving the necessary Qn values. Then the growth parts of the time histories would be found from Equation 9 up to about the inflection points of these S-shaped curves, and from Equation 10 thereafter. Since an exact solution is available for the first vessel, the approximations for ai are perhaps superfluous, but they may be used as a check against blunders and are included in the tabulation for completeness.
=
...
constant
= a.
Qn-1
n = 1,2,
+ an =
n = 1 , 2,
an-1 a0
=
. ..
(13)
constant
has the solutions
- e-T) az = a o [ l - (I + T)e-Tl a3 = a o [ l - ( 1 + T +T).-.] a1
=
ao(l
These solutions have the fault of tending to a0 for all tanks-which clearly does not happen when a reaction is taking place. T h e first suggestion was to replace the a0 in each of the Equations 14 by the appropriate equilibrium value U ~ E . Since the assumption of no reaction means t h a t a maximum amount of reactant A is entering each vessel, the approximations found must give values for an which are always too laree even when adjusted to the proper equilibrium concentration. The set of approximations that arises when the quadratic term is small is more complicated, but may be found by rather tedious
INDUSTRIAL AND ENGINEERING CHEMISTRY
Vol. 47, No. 4
ENGINEERING. DESIGN. AND PROCESS DEVELOPMENT VALUES USED IN ONE CYCLE
integration by parts. Equation I 1 is written with a new variable, S,which is equal t o B T , and the small but mathematically embarrassing terms are transposed t o the right. Thus,
F
COMPUTATION
d Q. dt
which is multiplied by e' and integrated
A fair approximation to a,(S) is made by ignoring the quadratic term in the integral, and then a much better approximation may be obtained by substituting this first approximation into the integral for a,. T o illustrate, the first reactor ( n = 1) is taken, where the first approximation is simply e"&
a0 = - (1
B
o
I
2
3
4
0 - POINTS KNOWN AT
START OF CYCLE
X = TENTATIVE POINT COMPUTED IN MIDDLE OF CYCLE
Q"
- e+)
t
O = F I N A L POINT OF CYCLE
which, on substitution into the Ka: term in Equation 16 gives a second approximation
al(S)= a B0 ( l
-
e-8)
-
(2) [I
2Se-S
- e--2s]
(17)
where S = BT. Since Equation 17 is the first two terms of an infinite series which usually converges, i t has the proper growth shape, but it tends to the wrong limiting value because the rest of the terms are missing. Accordingly, the coefficients are adjusted so as to give the proper asymptotic values. When this is done, the equation becomes
Figure 2.
Sequence of numerical integration 2a0
=B
+ & coth ( & T / 2 )
the relationship can be written either as as given with the main body of the formulas. Substitution of Equation 17 into Equation 16 will yield a two-term first approximation for a2, which will yield a five-term second approximation. This process rapidly becomes unwieldly, and terms which involve the initial concentration, ao, to the first and second powers have been retained on the grounds t h a t if a0 is small the others are unimportant. However, the higher terms actually involve Kao, which may not be small. If Kao is one or larger, these solutions are poor approximations.' I n general, these approximations, when adjusted to have the correct limiting values, are good for the-early parts of the rise curves-Le., when t is small and the inflection point has not yet been reached. This technique of successive approximations has given the same solutions, Equation 9, that may be obtained by the formal reversion-of-series method for nonlinear differential equations. The other useful approximations, Equation 10, were suggested by the exact solution to the first reactor. If a l E is introduced into
Table 1.
a1
e-QT)
if B / & = 1
(20)
al =
alE
tanh(QT/Z)
if B / Q = 0
(21)
I n practice, B / &seems to lie between 1/2 and 1, so that Equation 20 is preferred. This suggested that the exponent, QT, be tried in earlier approximations, Equation 14, instead of BT. The formulas, Equation 10, resulted and have proved satisfactory for larger values of t (above the inflection point) when they were adjusted by replacing the multiplicative constant with the correct limiting value. Numerical integration
T o integrate Equation 3 numerically, the I B M card-programed calculator of the Forrestal Research Laboratory a t Princeton, N. J., was used. Milne's method in integration was
Agreement between Correct Answer and Approximations
Basic Parameters
a2
a2
aa
K
E
B
ao
B/Q
20
0.2
5
0.15
0.8
Correct value Error from Equation 9 Error from Equation 10
0.001982 4 -96
1
0
1
0.06
0.9
Correct value Error from Equation 9 Error from Equation 10
0.00608 -1 33
20
0
1
0.15
0.3 0.5
Correct value Error from Equation 9 Error from Equation 10
0.00358
1.0
= alE(1 -
or as
t = 0.24 0.000142 0 - 13
t = 0.60
t
= 0.24
0.000653 3 -38
0.03806 0 19
0.01544 217 100
a3
a3
a2
0.00529 1 -1
0.001046 -9 -5
t = 6.84
0.027408 140 270
-
0.04542 3 0
0.00310 28 -11
0.03721 31665 8
t = 0.60
0.00028 -2 4
8
0.00446 3 - 12
a2
t = 1.56
t = 3.00
0.001153 -2 99
94
a3
t = 0.60
0.042813 51 - 14
t = 3.00 0.02422 30945 87
t = 5.88
0.03731 863" 0
0.02490 1768" -1
"Approximations (Equation 9) are worthless for large t and these extreme values of the physical parameters; however, the coefficient of suitability is small. ~~~
April 1955
~
~
~~~
INDUSTRIAL AND ENGINEERING CHEMISTRY
~~~
~~
709
ENGINEERING. DESIGN. AND PROCESS DEVELOPMENT da are evaluated at dt regularly spaced intervals. Although the procedure is numerical i t may easily be described in geometric terms. A parabola is da fitted to the latest three values of and extrapolated to allow
used (IO).
I n this method, both a, and
_f
the area under it from to to t, to be obtained, as shown by the da, from dt t o to t 4 and hence is the increment to be added to a,(t,) to get an approximate value of an(t4). This approximate value of a, is then substituted in the differential equation, Equation 3, to give da, da a value for - a t t4. With this value of 2,a new and more dt dun accurate integration of - is made-this time from t P to t 4 , as dt shown by the crosshatched area in Figure 2. This incremental area, As, is to be added to a,(t2) to give an improved value for 4 1 4 ) . If the agreement between the two values of an a t ta is satisfactory (within approximately 20 units in the least significant figure that is to be kept a t each step), the values of an and dun - are considered to be known a t t 4 , and the process is advanced dt one step. Starting the Milne method requires four good consecutive values for amand three for its derivative, which may be obtained by special but analogous formulas given in the standard accounts of the method. The I B M card-programed calculator has sufficient internal da storage to keep the four values of an and the three of $ for each shaded area of Figure 2. This area, Ai, is the integral of
of two reactors (No. 2 and No. 3, since S o . 1 is computed exactly), as well as the current value of t d and its increment, At, each in its characteristic memory location. An instruction deck of cards was then prepared to carry through one cycle of the Milne method, and each function was finally replaced by its successor in the time scale. Thus iteration of the instruction deck causes the integration calculation to proceed. By watching the agreement of the two values for the a , generated a t each step of the computation, the operator can decide whether the time interval, At, is too small or satisfactory. If the agreement becomes inefficiently good, the operator doubles the interval (which is necessarily small in the beginning), repeating the doubling whenever it seems justified by the demands of accuracy. I n addition to the Milne integration, three tentative approximations were coded into the routine-none of which proved very satisfactory. The best one turned out to be
a?i(t) =
B
+
2a,-i(T) Q n - 1 ~0th(&n-,T/2)
which was suggested by the closed form for the first reactor. Together with these tentative approximations, a whole set of curves for the three vessels could be Milne-integrated in about an hour. Without these unsatisfactory approximations this time could have been halved. . Later approximations, Equations 9 and 10, were found to be preferable and were coded in a separate routine. A whole set of curves for the three vessels could be obtained in 6 minutes. This saving largely reflects the fact that the approximations need be evaluated only a t the widely spaced points of interest (about six per curve), whereas to get those same variables by numerical integration they must also be evaluated a t the many intermediate points needed to prevent the integration error from accumulating.
710
Accuracy of approximations and useful range of parameters are examined
I n general, the larger B and the smaller K and a,the better these approximations become. A rough “index of suitability” for the approximations is B / & , which must lie between zero and one. If B / & is greater than 0.5, then these approximations exhibit (roughly) no worse than 10% error and usually are much better. For values of B / & about 0.9 or greater, the approximations are usually better than 0.5%. (This observation is no mathematical theorem, but i t is probably useful.) The data of Table I illustrate the type of agreement observed between the correct answer and the two approximations. For several values of the basic parameters, Table I gives the correct ) a,(t) with the errors of the two approximations. values for ~ ( tand The errors from Equation 9 and the discrepancy of Equation 10 are both given in units of the last digit for the corresponding a,. The range of the parameters explored in this work are l l n I 3 0.05 5 a. I0.15 O l K l 2 0
5 2 E 2 . 1 but that the limits of the left are not critical for the accuracy of the approximations. Nomenclature = volumetric rate of solvent flow (assumed constant), mi./ min a,, = weight fraction of reactant A in nth reactor, grams/gram b,& = weight fraction of reactant B in nth reactor, grams/gram V = vessel volume occupied by solvent medium, ml. k = reaction velocity constant, min.-’ I = time, min. E = absolute excess of reactant B above reactant A = bo ao, grams/gram 6 = V / q = nominal holding time in tank, min.
q
K = B = T = S = Qn = R, =
.
k6
l+KE t/8
BT
iB* + 4 K a , ] ’ / *
&,T
Subscripts n refers to reactor number = 1, 2, 3 . . 0 refers t o hypothetical reactor numher zero (corresponds to feed to reactor 1 ) Literature cited
(1) Denbigh, K. G.. T r a m . Faraday Soc., 40, 352 (1944). (2) Ibid., 43, 648 (1947). (3) Eldridge, J. W., and Pilet, E. L., Chem. Ew. Progr., 46, 290 (1950). (4) Golle, H. A., J . Agr. Food Chem., 1, 12 (1953). (5) Jones, R. W., Chem. Eng. Progr., 47, 46 (1951). (6) LeClerc, V. R., Chem. Eng. Sci., 2, 213 (1953). (7) MacMullin, R. B., and Weber, JI., Jr., Chem. Met. Eng., 52, 101 (1945). (8) MacMullin, R. B., and Weber, M., Jr., Trans. Am. Inst. Chem. Engrs., 31, 409 (1935). (9) Mason, D. R., and Piret, E. L , INL). ENG.CHEM.,42,817 (1950). (10) Milne, W. E., “Numerical Calculus,” Princeton University Press, Princeton, N. J . , 1949. (11) Perry, R. H., and Pigford, R. L., IND.ENO.CHEM.,45, 1247 (1953). (12) Weber, A. P., Chem. Eng. PTogr , 49, 26 (1953). RECEIVED for review July 10, 1964.
INDUSTRIAL AND ENGINEERING CHEMISTRY
ACCEPTED Ootober 2, 1954.
Vol. 47, No. 4
