Scale-up criteria for stirred tank reactors

Scale-up Criteria for Stirred Tank Reactors z J. J. EVANGELISTA, STANLEY KATZ, and REUEL SHINNAR T h e City College, CUNY, N e w York, N e w Yark zyxwvutsrqpon zyxwvutsrqpon A method is derived for the design of stirred tank reactors for homogeneous reactions. A simple mixing model proposed previously by Curl (4) is used to compute the effects of finite mixing time on complex chemical reactions. It is also shown how the parameters of the model can be obtained by tracer experiments, or estimated theoretically by the assumption of isotropic turbulence. It i s shown that in many practical cases the assumption of ideal mixing i s a good approximation. However, the effects of imperfect mixing are more likely to be felt in a large reactor than in a pilot plant. Some quantitative examples are discussed. Methods are derived to compute the average outlet concentration for complex systems such as autothermic reactions, polymerization, crystallization, etc. In a number of recent publications the effect of mixing on homogenous chemical reaccions is discussed. In this paper an attempt is made to summarize this work into a design and scale up procedure for stirred tank reactors. The method as described is limited to reactants which completely mix though it can be extended to certain heterogeneous reaction systems. Now in any such scale up problem the first question that one has to answer is what are the dangers involved. If there is only a simple one path reaction, then the only question is how much is t,he conversion affected, and this can be relatively easily handled either by safety factors or by some of the estimation procedures described in the following. The more complex and challenging case is the one in which the product quality itself might be affected by the scale up, and in this case there is no way of compensating for our lack of precise knowledge by a larger reactor volume. This is true especially in complex reactions, where side reactions might be favored by local overconcentrations. Other systems very sensitive to local overconcentration are systems involving nucleation, such as crystallization and certain polymerization processes. The problem of scaling up a homogeneous reaction is therefore basically to evaluate the effect of mixing on the reaction. The first question that we have to ask ourselves is what changes when we scale up. One property we want to maintain during scale up is a similarity in the basic flow regime. It was shown in previous work (8, 9) that this criterion can be fulfilled by choosing the size of the pilot plant or bench reactor such that the Reynolds number is large (> 104 or preferably 105). This ensures that the average velocity distribution is a function of space coordinates only and fairly independent of Reynolds number, resulting in a similarity in the overall velocity distribution between geometrically similar reactors. The second criterion commonly used and explained elsewhere (8, 9) is that the energy input per unit volume should be constant, as this gives a similarity of the turbulent flow regime in the high wave number range of the turbulent velocity spectrum. This is important for heterogeneous systems, as the velocity field that a single particle sees around its periphery remains constant during scale up. It was pointed out that the time scale of mixing changes during scale up at constant energy input per unit volume. To keep this time constant during scale up of geometrically similar vessels (with large Reynolds numbers) one would have to keep the agitator speed constant. For homogeneous reactions the overall mixing time is really the important criterion, and as long as we can keep both the Reynolds number high and the revolutions per minute constant, then we could scale up with considerable confidence. This is however in most cases imwactical. At constant revolutions per minute the Reynolds number increases linearly with the characteristic length L and the energy input per unit volume increases proportionally to Lz, If we want to keep the Reynolds number above 104 in the small vessel we either have a relatively small scale up ratio or we end up with an unrealistic energy consumption in the large vessel. We therefore often have to live with the fact that during scale up the mixing time increases and try to estimate this effect in a quantitative way. These considerations apply not only for the design of continuous reactors but also to the design of batch reactors, if one of the reactants is added continuously. Zweitering ( 1 7 ) has shown that if the reaction is of nth order and the feed is premixed one can estimate bounds on the conversion. One bound ( a maximum for n > 1 and a minimum for n < 1) is the case of maximum segregation, where all particles are assumed to mix only with particles having the same age and waiting time, the other bound is the case of maximum mixedness, which in the case of a stirred reactor with a Poisson residence time distribution is the same as t,he ideally mixed reactor. Obviously for a first-order reaction the conversion is independent of mixing and depends only on residence time distribution. For more complex reactions these limiting cases do not result in any bounds on conversion, and other methods to obtain these bounds for more complex systems have been proposed ( 1 5 ) . While these methods are not rigorous they still give a fairly good estimate of the bounds. Such bounding is especially useful if the bounds are close together. This indicates that the system is insensitive to mixing and can be scaled up quite safely. However, in the cases where it matters most, namely, in complex reactions, nucleation, etc., the bounds are very far apart, and thus bounding methods based on residence time distribution alone are not very useful. Luckily most agitated reactors are, in their behavior, very close to ideally stirred tanks. Mixing times in most reactors are measured in seconds whereas residence times are normally measured in minutes or hours. Unless we deal with a system sensitive to mixing or with very large reactors, it is very hard to measure any deviations from complete mixing [see for example (IS)].We therefore normally deal with the case of a reactor which at least in the pilot plant is very close to ideal mixing and our main problem is to estimate 1. How large is the deviation from ideal mixing likely to become in the large scale reactor? 2. How sensitive is the process to small deviations from complete mixing? A quantitative approach to these two problems is outlined in the following sections. zyxw zyxwvu zyxwv zyxwvutsrq Vol. 15, No. 6 THE MODEL The mixing model for a stirred tank reactor used here AlChE Journal Page 843 zyxwvuts zyxwvu is the one used by Curl ( 4 ) , Spielman and Levenspiel (11),Kattan and Adler ( 1 8 ) , and others in this and related contexts. Briefly, it regards the reacting mass as made up of a large number of equally-sized parcels of material (particles), which from time to time undergo independent pair collisions, equalize concentrations, and then separate. Between collisions, each parcel of fluid behaves like a little batch reactor. Fresh material is fed at a constant rate, and the withdrawal takes a representative cut of the contents of the vessel. One speaks of these particles as though they had a definite material identity. This is reasonable for dispersed phase systems but here where we are dealing with homogenous phase systems, they might perhaps be better regarded for the present purpose as primitive representations of turbulent eddies. If we consider, for concreteness, a single reaction where the concentration c of reagent behaves in batch according to x’ + x” .) (-y’- y” 2 y ) dx‘dy”d’dy’’ zyxwvutsrqp zyxwvutsr zyxwvuts zyxwvutsr zyxwvutsr zyxwvu zyxwvutsr then we may describe the contents of the reactor at any time t by the concentration distribution p ( c , t ) of particles at that time. Technically speaking, p is a probability density in c, with p(c, t ) d c giving the proportion of S, particles having concentration between a and b at time t. The distribution p, according to what has been said, satisfies the integro-differential equation ( 4 ) Calculations of reactor performance according to Equations (2) and (4) for some common one and two reaction schemes will be found below. The basic equation in p may be applied as well in the study of certain polymerization reactions, where we are concerned, to describe the molecular weight distribution. In a homogeneous radical polymerization, for example, we are concerned to describe the concentration of initiator (catalyst) i ( t ) , of monomer m ( t ) , of growing radicals dr, + ( r , t ) d r , having molecular weight between r and r and of terminated polymer having molecular weight between T and T + dr, + ( r , t ) d r . Assuming termination by combination, the batch performance of such a polymerization system may be described following reference 19 by + a$(r, t ) ar +Gm- at = d i(t)d(r) - D +(r, t ) ,f+dr =-Bi(t) where po(c, t ) is the concentration distribution for the feed, l/a the nominal residence time of material in the tank, and ,f3 a measure of the agglomerative mixing intensity. The mean concentration of reagent in the vessel and outlet is simply the first moment of p , and this is the working measure of overall reactor performance for a given kinetic system. Information about the random fluctuations in a turbulent mixing system is given by the higher moments. It should be noted that the concentration variable c in Equation ( 2 ) may be scaled in any one of a number of convenient ways, provided only that the basic property of averaging on agglomeration is preserved. Thus, we may write (2) with yield, conversion, extent of reaction, and so on, in place of c. Further, when there is more than one reaction going on, so that the kinetics are described by several simultaneous equations in place of Equation (I), we find natural generalizations of ( 2 ) . Thus, with two independent reactions whose progress variables (say) follow the kinetic scheme = - G m ( t ) S#Jdr Here, the molecular weight is reckoned in monomer units, v is the number of radicals formed by each molecule of decomposing initiator, and B, G, D are rate constants for initiation, propagation, and termination, respectively. Again, introducing the moments of the size distributions we find, (19) a set of batch kinetic equations in these moments and the concentrations d, m: dxo = uBi - DxO‘ dt dxi = Gmxodt DX&~ = 2Gmxl - D x g 2 dt d3c2 -= dyo dt D 2 -= dy2 D (xoxz xo2 (3) we have dt Page 844 AlChE Journal + x12) November, 1969 zyx zyz zyxwvutsrqp di -=-Bi dt drn dt -=- INTERPRETATION OF TRACER EXPERIMENTS zyx zyxwvuts Gntx, These serve as an eight dimensional version of ( 3 ) , and we have corresponding to them an eight dimensional version of (4) in the distribution p ( ~~,~~,~~,y~,y~,y2,$rn,t,i,m,t). From the distribution p , we may recover, say, the weight average molecular weight of radical together with polymer as Entirely similar considerations may be applied to the study of such particulate systems as crystallization reaction. Now the methods of solution of Equations ( 2 ) and ( 4 ) that are developed below apply regardless of the number of independent reaction variables, provided only that the kinetic expressions may be taken as polynomials in these reaction variables, They may accordingly furnish a useful guide to how the mixing interacts with strong nonlinearities in the kinetics, as in the crystallization nucleation rate, or with an independently varied critical feed stream, as in the initiator feed to polymerization reactors. We do not however present any results along these lines here, but reserve these studies for later report. Solutions to E uation (2) for several simple reactions have been publis ed ( 4 , 7, 11 ) In the literature (11) a Monte Carlo method was used and ( 7 ) a direct numerical method was employed. A method will be described in this paper which allows one to obtain the moments of p ( c ) in a simpler way. The authors want in no way to imply that the above model represents the real mixing processes in a turbulent reactor. These are far more complex and defy as yet an analytical description. The above model has, however, one important similarity to turbulent mixing. As will be shown below a concentration disturbance (in the absence of reaction) as computed from this model shows the same sort of decay as in isotropic turbulence, the variance (or the second moment) of the concentration in fluctuations in both cases decreasing exponentially with time. We furthermore do not claim that it is possible to predict a conversion accurately if B/a is not very large. Our claim that this simplified model is useful in design is based o n the fact that the behavior of the conversion with respect to B/a reaches an asymptotic value for high values of p/a. As long as B/a during scale up stays lar e enough so that we remain in a region where our resu ts are insensitive to p/a, we have good justification to assume that a real mixing process will also be insensitive to mixing in the same range of mixing times. In the design of real stirred tank reactors we mostly deal with small deviations from complete mixing. To say that the results are insensitive to the value of B/a provided B/a is large enough is the same as assuming the vessel is ideally mixed. The method described above allows one to estimate the mixing intensity required to approach ideal mixing. Now even if the vessel is not ideally mixed, and the results depend on p/a, as long as the deviations from ideal mixing are small, the above method should estimate these effects fairly accurately. The question that we still have to answer is how can we relate p/a to an experimentally measurable quantity, or how can we estimate it in the absence of suitable measurements. x We are concerned here to attach am empirical interpretation to the mixing intensity parameter B of Equation (2) and its higher-dimensional analogues. This interpretation is to be found in the analysis of tracer mixing experiments, and we accordingly begin by writing ( 2 ) , without the kinetic term, in the form We may note at once that the information we want is contained in the statistical fluctuations of concentration about its mean. As far as mean values go, the mixing system (without reaction) behaves like a completely mixed vessel with input-output time constant (Y. Indeed, if we denote the mean concentration by we find from ( 6 ) that equation p satisfies the ordinary differential (7) . independent of 8, where is the mean of the feed distribution po. And if we put a step of tracer into the feed of an initially tracer-free vessel, so that po(c, t ) = 6 ( c - c o ) ; t > 0 (8) and do) = 0 (9) z zyxwvutsrq zyxwvut zyxwvuts Y Vol. 15, No. 6 we find from ( 7 ) that the mean tracer concentration in the vessel and its outlet is p(t) =~ ~ ( 1 - e ~ “ ~ ) (10) the familiar stirred tank exponential response. T.he first measure of the concentration fluctuation is given by its variance V2(t) = J { c - P ( t H 2 p ( c , t ) d c and we find from Equations (6) and ( 7 ) that it satisfies the differential equation is the variance of the feed distribution. Again, with a step of tracer (8) in the feed of a vessel that is initially tracerfree, so that we have, besides, ( 9 ) , the initial condition AlChE Journal $(o) = 0 (12) Page 845 zyxwvutsrqp zyxwvutsrqpo zyx zyxwvut zyxwvu zyxwvut we find by solving (11) and consulting (10) that the variance of tracer concentration in the vessel and its outlet is given by g-at u"t) =a zyx zyxwvut - e-Pt co2 e - a t I3-a (13) With careful monitoring of an output line, the expression for the variance in (13) may serve as a guide to the estimation of /3. This variance rises from 0 to a maximum at 1 8+a t=In 8-. 2a - before decaying again to 0, and, with CY known, the bare location of this maximum will give some estimate of 8, even without detailed knowledge of the variance history. The step input experiments described above furnish perhaps the most convenient means for an empirical determination of p, although the underlying Equation (6) can readily be made to yield results appropriate to quite different experimental situations. Indeed, the most direct way of visualizing the effect of /3 is via a batch experiment with no input-output at all. If, in ( 6 ) , we drop the inputoutput terms, we are left with ( 4 ) S ( T - c ) and it is hoped that in the near future more measurements will be available, enabling one to correlate /3 with agitator design. I t may be noted finally that each value of the mixing intensity /3 corresponds to a definite value of Zweiterings degree of segregation (17) in the reactor. Indeed, if we denote the degree of segregation by s, we have simply dc'dc"-p(c,t) } a+l3 with vanishing @ giving complete segregation, and infinite 8, complete mixing. Thus, an estimate- of the degree of segregation, obtained perhaps by studying a particular reaction system (not first order), can be translated directly into an estimate of 8. (14) APPLICATION OF TURBULENCE THEORY The work of Batchelor (1) and Corrsin (3) and others on the application of the theory of turbulence to the scale up of stirred tank reactors leads directly to an estimate of the mixing intensity B in terms of the large-scale properties of the turbulence. Accordingly, we recapitulate the basic results here, and discuss their bearing on the practical estimation of 8. The analysis is all for a homogenous isotropic turbulence. The development begins with the standard partial differential equation for Fick's law diffusion superimposed on turbulent convection in an isolated system, from which one argues in the usual way that dc2 -=- dt dt and that the variance decays according to Thus, in a pure batch experiment, where a quantity of tracer is injected initially into some portion of the vessel, the resulting inhomogeneity in the distribution of tracer in the vessel decays, as far as variance goes, according to that is, simply with the time constant Js. According to this direct interpretation, /3 can be related to the characteristics of the turbulent flow in the vessel, and we shall review below tshese relations and their implication for the translation of j3 from one scale of operation to another. In order to get an experimental measurement of /3 we can therefore perform three types of experiments. One is to introduce an amount of tracer and measure the decay of the concentration fluctuations. The second is to introduce a step in tracer concentration and perform the same measurement. The third is in the case of multiple feeds to introduce a tracer in one of the feed streams only and measure the variance of the concentration fluctuations over the vessel at steady state. (The latter can also readily be computed.) In each case we need a method which is able to measure concentration fluctuations on a very small scale. Two methods have been recently developed which should make an important contribution to our understanding of the mixing process in stirred tanks. One developed by Kim and Wilhelm ( 5 ) is based on the scattering of light due to gradients in the refractive index. The second by Brodkey (2) uses fiber optics and colored tracers. As yet there have been very few actual studies 2 0 (grad c ) 2 Here D is the molecular diffusivity, c is a concentration fluctuation, and the overbar denotes an average so that c2 is simply the variance uz of concentration in a batch experiment. The mean square gradient in concentration fluctuation is related back to 02 by introducing a microscale 1, for the turbulent mixing, according to which the mean square components of the concentration gradient all have the common value d ( t ) = a(0)e-Pt Page 846 a zyx and we may see from this equation that the mean tracer concentration is constant -dP =o s=- 2 so that and dc2 -=-- dt 120 - L2 C2 This represents a simple exponential decay in variance, and consulting (15), we see that Information about the turbulent velocity field itself is however more naturally expressed in terms of the Taylor (dissipation) microscale l d , and of a corresponding Reynolds number R. The characteristic velocity appearing in this Reynolds number is the root mean square velocity fluctuation u. For an isotropic turbulence, this root mean square fluctuation is the same in each coordinate direction. AlChE Journal u5 = uy = uz = ib fl November, 1969 zyxwvutsrqpo zyxwvutsr zyxw zyxwvut zyxwvut zyxw and the appropriate Reynolds number is accordingly defined as R=-- 1 * * and, consulting (16) kdu (17) where v is the (molecular) kinematic viscosity of the working fluid, The two microscales are related (approximately) by the fact that this relation being valid for R large, but D / v not too large. We may note, as a guide, that difEusivities in water are of the order of to sq. cm./sec., while the kinematic viscosity v is 10-2 sq. cm./sec. In most practical liquid reaction problems D / v is less than or equal to unity. [Should D / v be large the literature (1 ) shows how the following analysis might be modified.] Now the Taylor microscale $I may be related to a characteristic linear dimension L of the mixing system as a whole (the paddle diameter, say, or the vessel diameter itself) by taking -l d= - A L R or, applying the definition (17) of the Reynolds number The first inference commonly drawn from ( 2 2 ) is that in order to keep the same reactor mixing properties on scale up, one should maintain similarity, and also keep ~ In geometrically similar vessels, the ratio E / L constant. as e is proportional to N3L2, this means keeping the agitator speed N constant during scale up. Also, since c is the power input per unit mass, keeping e/L2 constant amounts to making the actual power input proportional to L5, a demand which often results in uncommonly high power requirements. More realistically, we may use ( 2 2 ) to give a numerical estimate of B for a specified reactor size L and power input per unit mass t, and then use the methods developed below to estimate the effect of this B on the reaction system in hand. To do this requires of course a knowledge of the constant A f , but since this is scale-free, it can be estimated economically in small scale by means of the tracer experiments discussed above. The dependence of on L and L shown in ( 2 2 ) might of course also have been obtained by a straightforward dimensional argument. But the analysis based on turbulence theory leads also to some idea of the actual numerical magnitudes involved that may serve as a useful guide to practice in the absence of the tracer experiments suggested above. From work in wind tunnels and in pipe flows, Corrsin cites the (approximate) values A FJ 2 0. ., ~ 1/2 FJ so that ( 2 2 ) gives Here A is an empirical parameter, independent of the scale of the mixing system, but varying from one mixing configuration to another. One would expect A to depend, for example, on the physical properties of the working fluid, on the shape of the stirrer in a mixed vessel, and on the ratio of stirrer diameter to vessel diameter, but to have the same value for geometrically similar mixing systems working on the same fluid. The Taylor microscale l d may also be related to the power input to the turbulent fluid. Using a conservative numerical factor, one takes the energy dissipation per unit mass of fluid to be 10/3 v u2/I?d2, and, with an empirical efficiency factor 7, one may identify this as the actual power input t to the mixing system per unit mass of working fluid 10 7pu2 E = (20) 3 262 1 In a typical commercial application, with L = 100 cm., E = 2 x 104 sq.~m./sec.~ ( 1 h.p./100 gal. of water), Equation ( 2 3 ) gives /3 .63 Vsec. It may be noted that with these same numerical values, ( 2 1 ) gives for the root mean square velocity fluctuation, u 340 cm./sec., a plausible figure for the circulation velocity in the vessel (which in this context we identify with the turbulent velocity fluctuation). In agitated vessels both A and 7 should strongly depend on agitator design and one would need more accurate studies before one could really predict A and 7 on the basis of the geometrical design of the vessel. Tshere are several studies available on the affect of agitator design on overall mixing time. Van De Vusse ( 1 3 ) measured mixing time by observing the disappearance of large scale refractive index gradients via a Schlieren method. While this is not the same as the tracer experiment described in this section, we can still get an estimate of B, by writing - FJ zyxwvutsrq zyx zyxwvutsrq -- The efficiency 7 measures that fraction of the power input to the system which goes directly to the turbulent velocity field (and is only later dissipated as heat), as distinct from that fraction which goes directly to heat. One would expect 7 also to be independent of the scale of the mixing system, but varying from one mixing configuration to another. The basic physical relations in question are now contained in (17), (19), ( 2 0 ) , and these may be solved directly for u, &, 1,. One finds u= 0.27A2 (-> 1,= (80.0 A2q) Vol. 15, No. 6 B=- m 7m where the constant of proportionally m depends on the sensitivity of the method used to determine the time of complete mixing ( m should vary from 1 to 5 ) . Choosing m = 1 gives a real lower bound on /3, and this is exactly what is needed in our case. For a baffled completely stirred turbo or paddle mixer, Van De Vusse recommends a scale up equation derived by dimensional arguments which is exactly equivalent to Equation ( 2 2 ) . If one compares values of #? estimated from Van De Vusse's work to Equation ( 2 3 ) one finds that the values from ( 2 3 ) are too high. However, one should remember that Van De Vusse measured disappearance of a dye, and that by choosing B = 1/Tm7 we made a very conservative estimate of 8. In view of this zyxw (&)'I3 1 D3L2 (7 J AlChE Journal Page 847 zyxw zyxwvutsr zyxwvutsrqpo zyxwvut zyxwvutsrqpon zyxwvuts zyxwvutsrq zyxwvuts zyxwvuts the experimental results are really in surprisingly good agreement with our a prwri estimate. Van De Vusse’s data correlate well with (3) 1/3 B-0.1 As said before the constant in this correlation strongly depends on agitator design and should preferably be either measured or in the absence of more accurate data estimated from the literature ( 1 3 ) . In case a high value of 6 is desirable strong emphasis should be given to correct agitator design and correct placement of the feed pipes. Multiple agitators and multiple feed points can increase the value of 6 considerably. CALCULATION OF REACTOR PERFORMANCE We are concerned here to establish a systematic way of solving Equation (2) and its higher-dimensional analogues such as Equation (4) under conditions of steady state reactor operation. These equations seem in general to be very difhcult to solve. and since our interest commonly centers on high mixing rates, we proceed by expansion in powers of a//?.What results is a succession of linear integral equations in the successive contributions to the overall p , and while these don’t seem to be easily solvable either, they do lead to sets of linear algebraic equations in the moments of the successive contributions to p . The procedure closes, that is, it gives at each level of approximation a self-contained set of equations in a finite number of leading moments, provided the rate expressions ( 1), (3) and so on are polynomials. The resulting linear equations in the moments can be solved by standard means. The overall average reactor performance is found finally by adding up the contributions to the first moments of p . We begin with the one dimensional Equation ( 2 ) . We drop the time dependence, introduce a neutral variable x to stand for concentration, yield or whatever, and divide through by the input-output time constant to find for the distribution p If f is a polynomial of degree N , then (26) furnishes an overall relation among the first N moments of p . In particular, for a first order reaction with rate constant k, where x is the concentration of reagent so that f ( x ) = -kx/a, we find for the mean reagent concentration in vessel and outlet the standard result 1+- k a! (independent of the mixing intensity f i ) , where p1;o is the mean reagent concentration in the feed. In general, however, p (but not po) and its moments depend on the mixing intensity, and we proceed by expanding (24), (26), (27) in powers of l/X. We remind ourselves that in almost all practical applications A is huge as compared to unity, and we therefore expect such an expansion to converge fairly rapidly. We write formally and m , . - so that = pn(j) Jxn p(*)(x)dx (30) Bringing (28) to (24) and equating coefficients in gives for p ( 0 ) p‘0’ (x) - JJp‘O’ ( x ’ ) p‘”’ (x”) (Y for dl) x‘ + xn x 2 + A S J p ( x ’ ) p ( x ” ) a(?- x) dx‘dx” (24) Here po ( x ) represents the feed distribution, commonly 6 ( x - xo) for a suitable xo, f ( x ) is the rate expression normalized on a, and and for the higher p ( j ) zyxwvutsr is a dimensionless mixing intensity. The first result we want from (24) is an overall average material balance. If we denote the moments of p by fi = and those of po by CL~;O= s s - pu-1) (x) x n p ( x ) dx , xn po(x) dx ~=Pl;o-cL (34) and in the higher p ( j ) (26) We note also here that p ( x ) is to be a properly normalized probability density, so that Page 848 d -{ f ( x ) p ( + I ) ( x ) } ; i = 2 , 3 . . .. .. (33) dx Similarly, bringing (28) and (29) to (26) and equating coefficients in A gives in p(O) we find this material balance by multiplying (24) through by x and integrating -Jf(x)p(x) ) dx‘dx” pl(j)- Jf(x)p(j)(x)dx = 0; i = 1 , 2 , . ... (35) Finally from (27), we have that AlChE Journal November, 1969 while zyxwvutsrqpon zyxwvu zyxwvutsrqp zyxwvuts zyxwvu zyxwvutsrq zy zyxwvu po'o' = 1 po(j) = 0; j = (36) 1,2, . . . . (37) Now the equation (31) in p(O) simply represents the behavior of the perfectly mixed reactor (infinite #), and its solution is p'O)(x) = 8 [ x - p 1 ' 0 ' ] (38) which satisfies the normalization (36). The mean value p1(O) is determined by the overall material balance (34) - f [ p l ' O ' l = Pl;O pl'o' (39) With p ( 0 ) in hand, we may attack the equation (32) in the first correction p") by taking partial moments (30), that is, by multiplying through by successive powers of x and integrating. We find, taking due account of (37), linear algebraic equations in the p n C 1 )of the form moments can be found at will, this causes no special difficulty. Finally, having evaluated the partial moments for as many i as desired, we may see explicitly the effect of the mixing intensity on the overall performance of the reactor from (29). p1 = p1'0' + 1x -p1(1) + -1 A2 p1'2' + . .. (43) This completes our solution of the one-dimensional equation (24). Entirely similar considerations apply to the higher dimensional versions of ( 2 ) such as ( 4 ) . For the working equation, we have in place of (24) . - -pn;o + n[pl(o)]n-l - f [ . ~ 1 ' ~ )n ] ;= 2,3, . . . . (40) there being no information in the moment equations for n = 1. This loss of information is however made up by the overall material balance (35) for j = 1. - s pl(l) S . zyxw (41) and taking it together with the Equations (40) for n = 2 , 3 , . . . ,N, we find a set of N linear algebraic equations in ~ l ( ~. ).,, p N ( l ) which can be solved by standard numerical means. With these first N moments in hand, as many further moments as desired can be generated by taking (40) for n = N 1, N+ 2, etc., in turn. The same procedure applies to the higher p ( j ) of (33). Taking moments, we find a set of algebraic equations very much like (40), indeed having the same left hand sides p ~ 1 ~ 1 x +2 ~ X2 R ~ R ( X) dx .. nR,o= Jx1nlxn2. . .XnR p o ( x ) c l j , There are now R overall material balances in place of (26) - Jfl(x)p(x)dx . + . . nR = pn1nZ with pn1nz f ( x ) p ( l ) ( x ) d x= 0 If f is a polynomial of degree N, Equation (41) is another independent relation among the first N moments of pcl), . where x stands for (XI, xz, . . , x R ) , dx for the volume element dxl, dxz, . . . , d X R , and the 6-symbol for the appropriate product of &functions. The moments of p are now multidimensional, and we denote them - = p1o Jf I t ( x ) p ( x ) d x = Po0 . . O;O - P I O . . . o . . .1;o-poo.. 1 (45) and these determine entirely the average behavior of the system for first order kinetics. The single normalization condition (27) stands. Tqheexpansion (28) applied to (44) leads to (31) for , A for p ( l ) , and to ... i = 2,3,. . . . n = 2,3,. (42) the right hand sides for each i involving moments of the earlier p ( j ) . Again, there is no information in the moment equations for n = 1, but i f f is a polynomial of degree N , we find for each i a closed set of equations in p l ( j ) , pz(j), . . , p ~ ( 5by ) taking Equation (42) for n = 2, 3 , . . . , N together with (35). As before, once the first N moments of p'i) have been found, as many higher moments as desired can be generated by taking (42) for n = N I, N 2 etc. It should be noted that the integral term on the right hand side of (42) leads to a certain cascading in the number of earlier moments required as we go to successive values of i. Thus, with f of degree N, we need: 2N - 1 moments of p ( l ) to find N moments of p ( 2 ) ; 3N - 2 moments of p ( 1 ) to find (2N - 1) moments of p ( 2 ) to find N moments of pC3); etc. But since higher . + Vol. 15, No. 6 + for the higher p ( j ) . The moment expansions m - with (j) PninZ . . . nR = S ... dx X ~ ~ ~ X Z ~ ~R Zn ~ p (( X) j ) (49) applied to the material balance (45) give R equations in P'O' AlChE Journal p1;o:. 0 - Jfl(x)p'"(x)dx = p 1 0 . . .o;o Page 849 zyxwvutsrqponmlk zyxwvutsrqpo zyxwvutsrqpo zyxwvutsrqponm zyxwv zyxwvuts zyxwvuts zyxwv zyxw S f ~ ( ~ ) p ( ~ ) ( x ) d x = ~ ~ ~ . (50) ..i;o ~o:P:.i- and R equations in each of the higher p(j) -J Pld!'. .o /LOO(!), .1 - =0 f1(x)p"'(x)dx fR(X)p'')(X)dx = 0 j = 1,2,. . . (51) Finally, the moilleiit expansions applied to the normalization ( 2 7 ) give (0) ... o = PO0 and (j) POO. . . 0 = 0 ; (52) 1 i= I,%... (53) The Equation (31) in p(O) is solved just as for one dimension. We have (0) ...0). = S(xl--/LlO p'O'(X) * (0) *8(xR-P00...1) (S4) satisfying the normalization ( 5 2 ) , with the R first moments determined by the R material balances (50). (0) PlO (0) . . . 0 - f l C P l 0 . . . 0, (0) (0) POO.. . 1 - f R [ b ' + l O . . . 0, . * . (0) * ,poo.. (0) 1 ? POO. Fig. 1. Second order reactions; 2A values. 3 - - - - - - - Asymptotic 28. (j) Pnin2.. nR- 2nl+ . . . + n ~ - l ml=O .ll = P l O . . .o;o . .11 = P O O . . . l;O (55) For the first correction p"), we take moments. We find from (46) (1) Pnln2. ' n -~2 n l + . . + n R - l ml=0 r=l 111 + n2 + . . . + j = 2,3,. . . nR >1 . (57) and for polynomial f r these equations can again be closed off by applying the overall balance (51) and the normalization condition (53). Tthere is again, as we take successive values of j , a cascading requirement for the number of moments of the earlier pj, but these can as before be generated in turn once each basic set of equations for the leading moments has been solved. Finally, having evaIu- Again, the leading term in the sum on the left hand side disappears by virtue of the normalization (53), and we find no information in the equations for n1 n2 . . . nR = 1. Assuming these f7 to be polynomials, each of degree L N in each argument, this deficiency is made up as in the 1 dimensional case by taking the Equations (56) with nl, n2, , . . , n R each running from 0 to N (but exn2 . . . n~ = 0, l ) , and adjoining to cluding nl them the R balance equations (51) for j = 1. What results is a self-contained set of N R R - 1 linear algebraic equations in the leading moments of p"). Even in complex cases, this is by no means a very large number for standard numerical methods. Thus, for the polymerization kinetics ( 5 ) , with N = 2 and R = 8, we would have only 23 equations. The procedure may as before be extended to the higher p'j) by taking moments in (47).We find + + + + + + + Page 850 + B + 2C. Separate feeds for A Fig. 2. Bimolecular reactions; A and B, pO(x, y) = 1/2[6 (x-2) = 1/2x0 = 1/2y0 = 1. - - AlChE Journal 6 (y) +6 ( 4 6 (y-2)1, - - - - - - Asymptotic values. CO November, 1969 zy zyxwv zy zyx not quoted here. In practical design problems X is almost invariably large as compared to unity, and the results for small values of A are therefore of little interest in our context. DISCUSSION AND EXAMPLES Fig. 3. Adiabatic reactions; x = concentration of product, yield = x averaged over contents of vessel. - - - - - - - - Asymptotic values. ated the partial moments for as many j as desired, we may see explicitly the effect of the mixing intensity on the overall reactor performance by assembling the R first moments from (48). 1 (0) P10.1 0 = PlO PO0 . l=POo 1 . . 0 +--lo.. h (1) .o +1 (2) 0 = p O . +-A1 ~ O(1). . . l + - 1P l o (0) ..l A2 (2) + , ,. o+ ... (58) This concludes our development of the solution to the reactor performance equations for high mixing intensity, and the series expansions and moment methods developed here are applied below to the calculation of reactor behavior for some common reaction systems. It might be noted that similar results can be developed for low mixing intensity by expanding the basic equations in powers of f3/a, that is, in powers of the A of ( 2 5 ) rather than 1/k. What results is a hierarchy of differential rather than integral equations in the successive contributions to p , the leading term being that corresponding to a completely segregated reactor. Moment methods are no longer appropriate, and indeed the successive (linear) differential equations can readily be solved: directly, for the onedimensional case; via characteristics, for higher dimension. The methods are quite straightforward, and the results are One can now summarize the preceding discussions into a design procedure for stirred tank reactors. The first and most important point is to determine whether deviations from complete mixing will have a detrimental effect. Quite often the opposite is true. Thus for example in a homogeneous premixed second-order reaction a degree of segregation different from zero would improve conversion. Several authors have therefore pointed out that assuming complete mixing will lead to an overestimate of the required volume. The authors do not agree to this as in our experience most agitated vessels are very close to complete mixing. It is not ractical to control the mixing so that it should be incomp ete, while maintaining good agitation heat transfer. If no agitation is necessary then we should a priori use a plug flow reactor in such a case. If mixing is necessary and it is further desirable to minimize the variance of the residence time distribution then this can be achieved by using a cascade of stirred tanks or similar devices. Therefore the main problem in designing or scaling up a stirred tank reactor is to decide whether small deviations from complete mixing have a detrimental effect or not. In Figure 1 the effect of A = 2f3/0r on conversion is plotted for some typical cases of a second-order reaction. We note that for values of A larger than 100 the conversion is already insensitive to variation of A. This is not the case lor A close to unity, but in stirred tank reactors such low values are quite uncommon. Consider again a practical example. For a 10 liter vessel with a standard turbine agitator and power input the low estimate of /3 according to the literature ( 1 3 ) would be 0.05 to 0.1 sec.-l. A value of X of 5 would mean a residence time of 50 sec. For a 10,000 liter vessel the lowest estimate of f3 under the same condition would be 0.01 to 0.02 and a value of A = 5 would correspond to 250 sec. residence time. As we stated already our main problem in the design and scale up of stirred tank reactors concerns reactions on which a too low value of X is detrimental. The simplest such case would be a reaction whose total order is less than unity. Such reactions are very scarce in liquid homogenous systems and again the effect of A is small.' Another fairly simple case is where two reactants enter the vessel in two separate feed pipes. In Figure 2 the conversion for a simple irreversible second-order reaction with separate feeds is plotted as a function of X. Again the effect is small for values of x over 100. Here however, in very large vessels the effect would lead to a significant reduction of A and should be evaluated quantitatively. A second case with similar features is an autothermic reaction, in which a cold feed is introduced into a hot reaction mixture. Again in this case incomplete mixing always reduces the conversion. In Figure 3 a numerical example of the effect of A on one typical autothermic reaction is given. In computing the effect of a on the reaction we made use of a method introduced by Spalding (lo), namely the fact that for most exothermic adiabatic reac- fl zy zyxwv zy Fig. 4. Parallel reactions; A + B; rate constant kl, 2A + 2C; rate constant kz. 5 = CB/CO-CA, selectivity. - - - - - - - Asymptotic values. - Vol. 15, No. 6 0 Several authors have discussed the case of a zero order reaction. The assumption of a zero order reaction is however a simplified description of the limiting reaction velocity of catalytic reaction at high pressures. These reactions are not zero order at high conversions, where the effect of A would be noticeable. Furthermore, our method (as well as Zwietering's) does not apply to heterogeneous catalytic reactions. AlChE Journal Page 851 zyxwvuts zyxwvutsrqpo zyxwvuts zyxwvutsrqpo zyxwvutsr tions the effect of conversion on reaction rate can with sufficient accuracy be described by -dx _ -A dt + (l-x)"(x E ) ~ (59) zyxw zyxwvutsrqp zyxwvutsrqponmlkj The constants n and m are determined by a best fit to the actual reaction rate dependence which normally contains an exponential term, such as dx -= dt A( 1- x) n c E I R T - x)ne-E/RTo(1 + ax) =A( 1 (60) Spalding has shown that for most exothermic reactions the mistake made in writing (59) instead of an Arrhenius relation is small. Equation (59) has the advantage that it can be treated by the moment method of the preceding section. A low value of leads always to a lower conversion and the effect is quite pronounced. However for values of a larger than 100 the effect of /3 again becomes small. These last cases are examples of a broad class of cases in which the only effect of a low value of a is a reduced conversion. We can always compensate for this by a larger residence time, and in fact we might want to compare the cost of intense mixing to the cost of an increased residence time. A quantitative evaluation, similar to the one in the examples given, will lead to a conservative design procedure as long as we use a conservative estimation procedure for A. A simple example of the second class of problems in which a low value of A might change the quality of the product is given in Figure 4, where the reaction is assumed to consist of two parallel reactions of different order [for an example see ( 1 4 )1. As the undesired reaction is of a higher order, the concentration of the reactant should be as low as possible everywhere. A low value of A will therefore favor the undesirable side reaction. If the side product is not separated, and the conversion is high the reduction of A due to scale up might have a serious effect. One notes here, that just increasing the residence time without changing the intensity of agitation will not be sufficient, as the side reaction will occur in the concentrated region near the inlet. There is a whole class of reactions of this type, and what complicates this case further is the fact that in many cases the exact mechanism and kinetics of all the possible side reactions is unknown. As another example of this type one might mention reactions in which the pH is controlled by acid addition and where high pH in the nonmixed region might cause side reactions. If the kinetics of all possible undesired side reactions is at least approximately known then the above methods can be used to estimate the minimum value of A necessary to guarantee the specifications. Sometimes it might turn out to be impractical to obtain a sufficiently high value of A in a very large tank and several parallel smaller reactors might provide simpler and cheaper solution. One might also want to investigate ways of increasing A by improving both the design of the agitator and the design of the feed distribution. It is unfortunate that we do not possess sufficiently accurate data for the effect of the feed distribution on agitator design. But it is apparent that using multiple agitators (single or multiple shaft) and multiple feed injection for each agitator, we can considerably increase A without increasing the energy consumption. Similar considerations apply to the case of systems with nucleation, in which the effect of A might be the most pronounced, as even with relatively very low concentration Page 852 fluctuations all the nucleation might occur in the small unmixed region. This will be discussed in detail in a future paper. One last and somewhat discouraging conclusion that we can draw from the above discussion is that direct empirical evaluation of the effect of agitator design on complex reactions is impossible in the pilot plant. If a reaction is sensitive to mixing time then this sensitivity is a strong function of /3 and therefore a reaction might be insensitive to agitator design in the pilot plant while it is very sensitive to correct agitator design in the large scale plant. This by the way is true of other systems sensitive to agitator design (9). ACKNOWLEDGMENT The work reported here was supported under AFOSR Grant No. 921-67. Some of the work reported here is a part of the research carried out by J. J. Evangelista in partial fulfillment of the requirements for the Ph.D. at The City University of New York. N OTATl ON A A = empirical parameter appearing in Equation (19) = rate constant in autothermic reaction expression a, b B = bounds on concentration = rate constant for initiation step in polymerization CA = outlet concentration of CB = outlet concentration of B = outlet concentration of C scheme cc A = concentration cf, Cf' = dummy variables co = concentration of A in feed = molecular diffusivity D = rate constant for termination step in polymerizaD C tion scheme f f, G = rate expression = rate expression, T = 1, 2 , . ,R = rate constant for propagation step in polymeriza- .. tion scheme i( t ) L = initiator concentration = reaction rate constant = characteristic length of the mixing system m m = microscale for turbulent mixing = constant of proportionality = exponent in approximate rate expression for auto- k & L, = Taylor (dissipation) microscale thermic reaction m ( t ) = monomer concentration N N n = agitator speed, rev./min. = degree of polynomial f = exponent in approximate rate expression for auto- n, = integer in moment definition; thermic reaction p po pj r = 1, 2, = concentration distribution ...R s = concentration distribution for the feed = jth term in expansion of p = Reynolds number = number of independent concentration variables = batch rate expression = linear dimension of crystal = molecular weight = degree of segregation s = batch rate expression t =time = root mean square of velocity fluctuation = x-component of u = y-component of u = z-component of u R R r r r u u, uy uZ AlChE Journal November, 1969 x x xj zyxwvutsrqpon zyxwvut zyxwvuts zyxwvutsrq zyxwvutsr zyxwvutsrqpo = Cartesian coordinate = progress variable = jth moment of size distribution of growing poly- mer chains x’, X” y y yj LITERATURE CITED = dummy variables = Cartesian coordinate = progress variable = jth moment of size distribution of terminated polymer y’, y” = dummy variables z = Cartesian coordinate Greek Letters a = inverse of nominal residence time /I =concentration distribution of growing polymer chains = concentration distribution of terminated polymer 1. Batchelor, G . K., J. Fluid Mech., 5, 113 (1959). 2. Brodkey, R. S., and J. 0. Nye, Reo. Sci. lnstr., 34, 1086 (1963). Lee, AIChE I., 10,187 (1964). 3. Corrsin: g:%d., 3, 329 (1957); 10,870 (1964). 4. Curl, R. L., AlChE I., 9, 175 ( 1963). 5. Kim, Y. G., and R. H. Wilhelm, Preprint 26B, Fifty-sixth Annual Meeting of the AIChE, Houston, Texas. 6. Rietema, K., Aduun. Chem. Eng., 5, (1964). 7. Shain, S., AlChE J., 12,806 (1966). 8. Shinnar, R., J. Fluid Mech., 10, Pt. 2, 259 (1961). , and J. M. Church, Ind. Eng. Chem., 52, 253 9. ( 1960); 53,479 (1961). 10. Spalding, D. B., Combustion Flume, 1,287,296 (1957). 11. Spielman, L. A,, and 0. Levenspiel, Chem. Eng. Sci., 20, 247 1965). 12. Tadmor, Z:, and A. Biesenberger, lnd. Eng. Chem. Fundamentals, 5, 336 (1966). 13. Van de Vusse, J. G., Chem. Eng. Sci., 4, 178,209 (1955). 14. lbid., 17,507,1962. 15. Weinstein, H., and R. J. Adler, Chem. Eng. Sci., 22, 65 (1967). 16. Worrell, G . R., and L. C. Eagleton, Can. I. Chem. Eng., 43,254 (Dec. 1964). 17. Zwietering, Th. N., Chem. Eng. Sci., 11, No. 1, 1 (1959). 18. Kattan, A., and R. J. Adler, AlChE J., 13,580 ( 1967). 19. Saidel, G. M., and S . Katz, ibid., 319 (1967). zyxwvu zyxwv = measure of agglomerative mixing intensity = actual power input to mixing system per unit mass of fluid 7 = empirical efficiency parameter A = mixing intensity parameter, = 2/3/a p = mean concentration = mean feed concentration ~ c n = nth moment of p f i ; o = nth moment of PO r-nf = nth moment of pj v = kinematic viscosity of workin fluid = number of radicals formed y each molecule of v decomposing initiator = variance of concentration distribution $ u.02 = variance of feed distribution T = mean residence time T~ = mixingtime c P; zyx Manuscript receitied October 11, 1967; retiision receitied June 26, 1968; paper accepted July 18, 1968. The Dynamic Modeling, Stability, and Control of a Continuous Stirred Tank Chemical Reactor W. FRED RAMIREZ and BRIAN A. TURNER University of Colorado, Boulder, Colorado A realistic CSTR model was developed and verified experimentally. The reaction studied was the exothermic, base sodium model was used to evaluate analysis, local linearizotion The effect of control valve hydroxide catalyzed decomposition of hydrogen peroxide. The the usefulness of the following stability analyses: steady state and Liapunov’s direct method through Krasovskii’s theorem. hysteresis on the system was olso investigated. The purpose of this investigation was to study both experhentally and theoretically the dynamics and control of a cantinuous stirred tank reactor, a C.S.T.R. There exist several analytical methods for determining stability in a C.S.T.R. (1 to 4, 7). These include steady analysis (71, hearization techniques (31, and Krasovskii’s Theorem with the Liapunov function it suggests ( 2 ) . Brian A. Turner is now with the Shell Oil Company, Martinez, Califomis. Vol. 15, No. 6 D Y N A M I C EQUATloNs The system considered was a continuous stirred reactor. The reaction was first order and homogeneously catalyzed. The reactant, hydrogen peroxide, and the cats$% sodium hydroxide, were fed into the reactor, while the product, water, and unreacted peroxide were continuously removed. The solution density and heat capacity varied with concentration, but their product is nearly constant and was considered constant in formulating the dynamic model. The exothermic heat of reaction was continuously removed through a cooling water belt which AlChE Journal Page 853