630
Ind. Eng. Chem. Res. 2003, 42, 630-635
GENERAL RESEARCH Analyses of Stability and Dynamic Patterns of a Continuous Crystallizer with a Size-Dependent Crystal Growth Rate Qiuxiang Yin,* Yanmei Song, and Jingkang Wang School of Chemical Engineering and Technology, Tianjin University, Tianjin 300072, People’s Republic of China
The effect of a size-dependent crystal growth rate on the stability and dynamics of an isothermal mixed suspension-mixed product removal (MSMPR) crystallizer is investigated. A linearized stability analysis has been performed to determine the stability boundaries in terms of dimensionless parameters that depend only on physical properties and operational variables. As in prior studies of crystallizer dynamics, the instability is always oscillatory and accompanied by a limit cycle. Numerical solutions of the actual nonlinear system equations are carried out to follow the cyclic behavior and to compute the amplitude of the limit cycle. that depend only on physical properties and operational parameters.
1. Introduction The stability and dynamics of continuous crystallizers have attracted considerable attention in the literature.1-7 Previous researchers have shown that instabilities in industrial crystallizers always manifest themselves in the form of limit cycle oscillations and that the stabilities of the crystallizers can be predicted using some relatively simple criteria. However, these criteria were derived by use of dimensionless groups that depend not only on the physical properties and control parameters in the problem but also on the steady-state values of the state variables. Jerauld et al.8 indicated that this would obscure the effect of the control parameters on the stability of the system. Therefore, they introduced three new dimensionless groups that depend only on physical properties and control parameters and derived a simple condition for the appearance of sustained oscillations in continuous crystallizers. From the same issues as considered by Jerauld et al.,8 similar dimensionless groups were recently used to study the multiplicities of the steady states in continuous isothermal crystallizers.9-17 Jerauld et al.8 used Volmer’s model and McCabe’s law to describe the kinetics of nucleation and crystal growth, respectively, but did not take into account the dependence of the crystal growth rate on the crystal size. It has been recognized that many systems exhibit anomalous growth with size dependency18-20 or growth-rate dispersion.21,22 Sherwin et al.1 reported that a size-dependent crystal growth rate has notable effects on the stability of continuous crystallizers. In the present work, the stability and dynamic behavior of an isothermal mixed suspension-mixed product removal (MSMPR) crystallizer with a size-dependent crystal growth rate is investigated. The stability boundaries are presented in terms of a group of dimensionless variables
2. Model Development An isothermal, unseeded MSMPR crystallizer of constant volume with a clear feed stream is considered. We assume that all new crystals are nucleated with infinitesimal sizes and that the breakage and agglomeration of crystals are negligible. Furthermore, the crystal growth rate is taken to be dependent on the particle size. Under these conditions, the population balance can be written as
∂n(r,t) ∂[G(r) n(r,t)] n(r,t) + )+ Bδ(r - 0) ∂t ∂r τ
(1)
where , the volume of liquid per unit volume of suspension, is related to n(r,t) and r by
∫0∞r3n(r,t) dr
) 1 - kv
(2)
Here, we consider the case of primary nucleation and represent the kinetic rate by Volmer’s model23
[
B ) kp exp -
kE (c/ce - 1)2
]
(3)
The size-dependence of the overall growth rate G is linear in the crystal size r and is modeled by the following expression1
G(r) ) kg(c - ce)(1 + ar)
(4)
A mass balance on solute leads to * To whom correspondence should be addressed. Tel.: (86)22-27405754. Fax: (86)22-27374971. E-mail: qxyin@ public.tpt.tj.cn.
d 1 [c + F(1 - )] ) [cf - c - F(1 - )] dt τ
10.1021/ie020493b CCC: $25.00 © 2003 American Chemical Society Published on Web 01/11/2003
(5)
Ind. Eng. Chem. Res., Vol. 42, No. 3, 2003 631
We define mi as the ith moment of n(r,t) through the following equation
mi )
∫0∞rin(r,t) dr
(6)
Substituting eqs 3 and 4 into eq 1 and recasting eqs 1 and 5 in terms of moments of n(r,t) lead to the infinite set of couple equations
[
]
m0 ke dm0 )+ kp(1 - kvm3) exp dt τ (c/ce - 1)2
assumption that we are investigating a first-order correction of the size-independent growth model, we estimate that β is small and cannot go beyond the range of (-0.3, 0.3). Among these dimensionless variables and parameters, Da, F, R, and β depend only on physical properties and operational parameters, and they can be varied by changing the inlet concentration and the residence time. Recasting eqs 7-9 into dimensionless form yields
dx0 ) -x0 + (1 - x3)Da exp(-F/y2) dθ
(10)
i ) 1, 2, ... (8)
dx1 ) -x1 + y(x0 + βx1) dθ
(11)
dc (cf - c) - 3kgkvτ(F - c)(c - ce)(m2 + am3) ) (9) dt τ(1 - kvm3)
dx2 ) -x2 + y(x1 + 2βx2) dθ
(12)
It is obvious that the infinite system of population balance equations for moments can be closed by the mass balance equation of solute at the equation for the third-order moment. Therefore, we will analyze the set of equations consisting of the mass balance and the first four moment equations. For ease of analysis, the following dimensionless variables and parameters are introduced:
dx3 ) -x3 + y(x2 + 3βx3) dθ
(13)
dy 1 - y - (R - y)y(x2 + 3βx3) ) dθ 1 - x3
(14)
mi dmi )+ ikg(c - ce)(mi-1 + ami) dt τ
(7)
Da ) 6kvkpkg3τ4(cf - ce)3 F ) kEce2/(cf - ce)2 R ) (F - ce)/(cf - ce) β ) akgτ(cf - ce)
In consequence, the stability and dynamic behavior of the crystallizer are governed by the dynamic model in eqs 10-14. 3. Linearized Stability Analysis At steady state, eqs 10-14 can be algebraically reduced to a single nonlinear equation in y as follows
(1 - ys)(1 - βys)(1 - 2βys)(1 - 3βys) ) (R - 1)Days3 exp(-F/ys2) (15)
θ ) t/τ x0 ) 6kvkg3τ3(cf - ce)3m0 x1 ) 6kvkg2τ2(cf - ce)2m1 x2 ) 3kvkgτ(cf - ce)m2 x3 ) kvm3 y ) (c - ce)/(cf - ce) Notice that the three dimensionless parameters, Da, F, and R, are identical to those of Jerauld et al.8 Da is a measure of the nucleation and growth that occurs in one residence time in the crystallizer, that is, as Da is increased, a greater mass of crystal is formed in one residence time; F can be interpreted as a measure of the nonlinearity of the primary nucleation rate and is also a prime determinant in the gross level of primary nucleation. The change in volume due to crystallization is measured by R; as R becomes larger, the volume change upon crystallization increases. The typical ranges of Da and R were estimated by Jerauld et al.8 as 0.110 000 and 10-100, respectively. We estimate that F is far smaller than 0.5 industrially, whereas it might be as high as 1 experimentally. β can be interpreted as a measure of the relative size increase in one residence time due to the size-dependent term ar. With the
Jerauld et al.8 has shown that, when β ) 0, eq 15 has a single root for all F, Da, and R; more generally, Lakatos14 indicated that this is also valid for all β. Therefore, the system under consideration has a unique steady state over the entire possible region of parameters. To assess the stability of the singular steady-state solution of the system of equations, we perform a general local stability analysis24 through the linearization of the model eqs 10-14 around the steady state and examine the eigenvalues of the Jacobian matrix. At steady state, the Jacobian matrix reduces to
[
J) -1 0 ys 0 0 0
0
-1 + βys 0 ys -1 + 2βys ys 0 (R - ys)ys 0 1 - x3s
-
x0s 1 - x3s
0 0 -1 + 3βys 3β(R - ys)ys 1 - x3s
2Fx0s ys3 x0s + βx1s x1s + 2βx2s x2s + 3βx3s (R - ys)(x2s + 3βx3s) -1 1 - x3s
]
(16)
The characteristic equation of the Jacobian determinant is
(-1)5 det(J - λI) ) (λ + 1)(λ4 + w1λ3 + w2λ2 + w3λ + w4) ) 0 (17)
632
Ind. Eng. Chem. Res., Vol. 42, No. 3, 2003
where
w1 ) 4 - 6βys + d
(18)
w2 ) 6 - 18βys + 11β2ys2 + (4 - 6βys)d
(19)
w3 ) 4 - 18βys + 22β2ys2 - 6β3ys3 + (6 - 18βys + 11β2ys2)d (20)
(
w4 ) (1 - 6βys + 11β2ys2 - 6β3ys3)
)
R - ys 2Fd + 2 + R-1 y s
(3 - 12βys + 11β2ys2)d (21) with
d)
(R - ys)(1 - ys) (R - 1)ys
>0
(22)
The above characteristic equation, although more general, is somewhat similar in structure to that obtained by Jerauld et al.8 for the case of size-independent crystal growth rate. The present result is valid for both sizeindependent and size-dependent crystal growth rates. According to the stability theory, to keep the system stable, the real parts of the eigenvalues of the characteristic eq 17 must be negative. Applying the classical Routh-Hurwitz criteria, we find that, for stability, the following conditions must hold:
Figure 1. Stability boundary as a plot of Da versus R for various values of F (β ) 0.15).
(i) All of the coefficients of the characteristic equation must be of the same sign, i.e., w1, w2, w3, and w4 should be positive. (ii) w1w2 - w3 > 0
(23)
(iii) w1w2w3 - w23 - w21w4 > 0
(24)
From the practical considerations, the values of w1, w2, w3, and w4 are bounded by the inequality constraints of 0 e ys e 1 and β < 0.3. It is easy to demonstrate that the first two conditions are always satisfied; therefore, stability can be determined from the third condition only. 4. Results and Discussion 4.1. Stability. In the previous section, we mentioned that the system of eqs 10-14 has a unique steady-state solution over the entire feasible region of parameters. This can also be verified by the fact that the eigenvalue of the characteristic eq 17 is never equal to zero, because w4 is never zero unless the value of β is greater than 1, which is practically impossible. Consequently, according to bifurcation theory,24 the real bifurcation of the steadystate solution of eqs 10-14 can never take place; thus, no change in the number of steady-state solutions can be expected. The stability boundary for the system can be obtained from the inequality in eq 24, which gives the conditions for the real part of the eigenvalues to be negative and the system to be stable. When the inequality is not satisfied, at least one of the eigenvalues is real and positive, or there is at least one pair of complex eigenvalues with a positive real part, and the system is unstable. Figures 1 and 2 show the stability boundaries in terms of a plot of Da versus R for different values of
Figure 2. Stability boundary as a plot of Da versus R for various values of β (F ) 0.5).
F and β. These stability boundaries are given only in terms of independent parameters and do not involve the steady-state variables. An interesting phenomenon, which can be seen from Figures 1 and 2, is that the curve separating stable and unstable regions plotted in R and Da forms a straight line in the logarithmic coordinate system. This was also shown by Jerauld et al.8 for the case of a size-independent growth rate (β ) 0). From the linearity of the plot and the slope of -1, it is possible to present the conditions for stability in a simpler manner; that is, we can produce a plot of a curve separating stable and unstable regions in F and the composite group RDa for different values of β, as shown in Figure 3. It is shown in Figure 3 that, for sufficiently small F or RDa, the system is always stable. It is also seen that a growth kinetic model with larger β raises the stability limits, leading to more possibilities for stable situations. Furthermore, increasing the residence time and keeping all else constant corresponds to increasing RDa at constant F (for example, from point A to B in Figure 3); it is clear from Figure 3 that increasing the residence time can cause instability. Changing the inlet concentration cf and keeping all else
Ind. Eng. Chem. Res., Vol. 42, No. 3, 2003 633
Figure 3. Stability boundary as a plot of RDa versus F for various values of β.
Figure 4. Projection of state-point trajectories behaves as a stable focus as Da changes from 800 to 500 (F ) 0.5, R ) 20, β ) 0.15).
constant corresponds to moving the operational point along a line RDaF ) constant (for example, from point C to D in Figure 3), where low inlet concentration corresponds to large F. From Figure 3, we see that with a decrease in cf, a stable system can become unstable. Therefore, an unstable system can be made stable simply by decreasing the residence time or increasing the inlet feed concentration, keeping all other parameters constant. 4.2. Hopf Bifurcation and Dynamic Behavior. We note that all of the coefficients of the characteristic polynomial are positive. From the Descartes rule of signs, the eigenvalues can never be positive and real. Therefore, the eigenvalues are a mixture of negative real, negative complex, and positive complex roots. Consequently, as shown in Figures 4 and 5, when the system is stable, it behaves as a stable focus, and when the system is unstable, it behaves as a stable limit cycle. Using the Hopf bifurcation theorem, Ferauld et al.8 showed that, in the case of a size-independent crystal growth rate, the instability always results in a Hopf bifurcation and hence is accompanied by a stable limit cycle. Because of algebraic complexity, the applicability
Figure 5. Projection of state-point trajectories behaves as a stable limit cycle as Da changes from 2500 to 1500 (F ) 0.5, R ) 20, β ) 0.15).
Figure 6. Bifurcation diagram of x3 versus Da for changing residence time τ (F ) 0.5, R ) 20, β ) 0.15).
of the Hopf bifurcation theorem could not be confirmed in this general case. However, solution of the dynamic eqs 10-14 indeed shows that the instability is always oscillatory and manifests itself as a stable limit cycle, as shown in Figure 5. Figure 6 illustrates the bifurcation diagram of x3, that is, the solid fraction (1 - ), versus Da for changing residence time τ (corresponding to moving the operational point along the line AB in Figure 3), showing the amplitude of the limit cycle around the unstable steady state represented by the dash-dotted line. With an increase in Da (corresponding to an increase in τ), the stability of the system is changed at the bifurcation point, and the amplitude of the oscillations increases as the system moves further away from the bifurcation point. The bifurcation diagram of x3 versus F for changing inlet concentration cf (corresponding to moving the operational point along the line CD in Figure 3) shows similar behavior, as depicted in Figure 7. In Figure 7, with an increase in F (corresponding to an increase in cf), an instability in conjunction with the appearance of a limit cycle appears at the bifurcation
634
Ind. Eng. Chem. Res., Vol. 42, No. 3, 2003 B ) nucleation rate, s-1‚m-3 c ) concentration, kg‚m-3 d ) dimensionless parameter defined by eq 22 Da ) dimensionless parameter ) 6kvkpkg3τ4(cf - ce)3 F ) dimensionless parameter ) kece2/(cf - ce)2 G(r) ) growth rate of crystals of radius r, m‚s-1 I ) identity matrix J ) Jacobian matrix kE ) constant in eq 3 kg ) growth rate constant, m4‚s-1‚kg-1 kp ) nucleation rate constant, s-1‚m-3 kv ) volume shape factor mi ) ith moment of n(r,t) n(r,t) ) population density at radius r ant time t, m-4 r ) crystal size, m t ) time, s w1, w2, w3, w4 ) coefficients of characteristic eq 17 xi ) ith dimensionless moment of n(r,t) y ) dimensionless concentration
Figure 7. Bifurcation diagram of x3 versus F for changing inlet concentration cf [R ) 25 × (10F)0.5, Da ) 4000 × (10F)-1.5].
point. With a further increase in F, the amplitude of the oscillations increases. It is well-accepted1,5,8 that the dominating source of this periodic behavior is the strongly nonlinear nature of the dependence of nucleation rate on supersaturation. Speaking in detail, a small increase in supersaturation will cause a temporary increase in nucleation rate, and then an excess amount of nuclei will grow, supplying such a large surface area for growth that the decreasing supersaturation causes the nucleation rate to decrease until the removal of crystals (and therefore surface area for growth) by convective outflow and the addition of more solute in the feed stream cause another increase in supersaturation. However, this effect occurs with a considerable time delay because the new nuclei have no appreciable surface area for a considerable time span. A kinetic model that predicts increasing growth rate with increasing size leads to a shorter time lag between the appearance of increased nuclei and increased surface area. Therefore, it can be expected that a model with a larger value of β should raise the stability limits of the system. This is consistent with the results displayed in Figure 2. 5. Conclusion We have developed a model for investigating the influence of a size-dependent crystal growth rate on the stability and dynamics of a continuous MSMPR crystallizer. Instabilities are possible over certain ranges of kinetic, physical, and operational parameters. The stability boundaries are defined in terms of a group of dimensionless variables that depend only on physical properties and operational parameters. The dynamic analysis indicates that stability is accompanied by a stable focus, whereas instability is accompanied by a stable limit cycle. In the regions close to the bifurcation point, the system exhibits extreme parametric sensitivity. Although the parameter ranges within which the instability might occur are narrow and the phenomenon of the oscillation is limited to some localized situations, a prior awareness of these regions is important in the design, start-up, and control of a crystallization process. Nomenclature a ) Constant in the size-dependent growth model equation (eq 4)
Greek Letters R ) dimensionless parameter ) (F - ce)/(cf - ce) β ) dimensionless parameter ) akgτ(cf - ce) δ ) Dirac delta function ) volume of liquid per unit volume of suspension λ ) eigenvalue of the characteristic matrix θ ) dimensionless time F ) crystal density, kg‚m-3 τ ) mean residence time, s Subscripts e ) saturation f ) feed s ) steady-state value
Literature Cited (1) Sherwin, W. B.; Shinnar, R.; Katz, S. Dynamic Behavior of the Well-Mixed Isothermal Crystallizer. AIChE J. 1967, 13, 1141. (2) Nyvlt, J.; Mullin, J. W. The Periodic Behavior of Continuous Crystallizers. Chem. Eng. Sci. 1970, 25, 131. (3) Lei, S. J.; Shinnar, R.; Katz, S. The Stability and Dynamic Behavior of a Continuous Crystallizer with a Fines Trap. AIChE J. 1971, 17, 1459. (4) Randolph, A. D.; Beer, G. L.; Keener, J. P. Stability of the Class II Classified Product Crystallizer with Fines Removal. AIChE J. 1973, 19, 1140. (5) Yu, K. M.; Douglas, J. M. Self-Generated Oscillations in Continuous Crystallizers, Part I: Analytical Prediction of the Oscillating Output. AIChE J. 1975, 21, 917. (6) Song, Y. H.; Douglas, J. M. Self-Generated Oscillations in Continuous Crystallizers, Part II: An Experimental Study of an Isothermal System. AIChE J. 1975, 21, 924. (7) Bhatia, S. K. Dynamics of Continuous Precipitation under Non-Stoichiometric Conditions. Chem. Eng. Sci. 1989, 44, 751. (8) Jerauld, G. R.; Vasatis, Y.; Doherty, M. F. Simple Conditions for the Appearance of Sustained Oscillations in Continuous Crystallizers. Chem. Eng. Sci. 1983, 38, 1765. (9) Tavare, N. S.; Garside, J. Multiplicity in Continuous MSMPR Crystallizers, Part I: Concentration Multiplicity in an Isothermal Crystallizer. AIChE J. 1985, 31, 1121. (10) Tavare, N. S.; Garside, J.; Akoglu, K. Multiplicity in Continuous MSMPR Crystallizers, Part II: Temperature Multiplicity in a Cooling Crystallizer. AIChE J. 1985, 31, 1128. (11) Tavare, N. S. Multiplicity in Continuous Crystallizers: Adiabatic Reactive Precipitation. Chem. Eng. Commun. 1989, 80, 135. (12) Padia, B. K.; Bhatia, S. K. Multiplicity and Stability Analysis of Agglomeration Controlled Precipitation. Chem. Eng. Commun. 1991, 104, 227. (13) Lakatos, B. G. Influence of Nucleation Mechanisms on the Multiplicity of Steady States in Isothermal CMSMPR Crystallizers. React. Kinet. Catal. Lett. 1992, 47, 305. (14) Lakatos, B. G. Uniqueness and Multiplicity in Isothermal MSMPR Crystallizers. AIChE J. 1996, 42, 285.
Ind. Eng. Chem. Res., Vol. 42, No. 3, 2003 635 (15) Yin, Q. X. Multiplicity in Continuous Adiabatic MSMPR Reactive Precipitators. Chin. J. Chem. Eng. 1998, 6, 138. (16) Yin, Q. X.; Wang, J. K.; Xu, Z.; Li, G. Z. Analysis of Concentration Multiplicity Patterns of Continuous Isothermal Mixed Suspension-Mixed Product Removal Reactive Precipitators. Ind Eng. Chem. Res. 2000, 39, 1437. (17) Yin, Q. X.; Wang, J. K.; Zhang, M. J.; Wang, Y. L. Influence of Nucleation Mechanisms on the Multiplicity Patterns of Agglomeration-Controlled Crystallization. Ind. Eng. Chem. Res. 2001, 40, 6221. (18) Abegg, C. F.; Stevens, J. D.; Larson, M. A. Crystal Size Distribution in Continuous Crystallizers when Growth Rate Is Size Dependent. AIChE J. 1968, 14, 188. (19) Mullin, J. W.; Jones, A. G. Programmed Cooling Crystallization of Potassium Sulphate Solutions. Chem. Eng. Sci. 1974, 29, 105. (20) Garside, J.; Jancic, S. J. Growth and Dissolution of Potash
Alum Crystals in the Subsieve Size Range. AIChE J. 1976, 22, 887. (21) Janse, A. H.; De Jong, E. J. The Occurrence of Growth Dispersion and Its Consequences. In Industrial Crystallization 1975; Mullin, J. W., Ed.; Plenum Press: New York, 1976. (22) Randolph, A. D.; White, E. T. Modelling Size Dispersion in the Prediction of Crystal Size Distribution. Chem. Eng. Sci. 1977, 32, 1067. (23) Nyvlt, J.; Sohnel, O.; Matuchova, M.; Broul, M. Kinetics of Industrial Crystallization; Elsevier: Amsterdam, 1985. (24) Kubicek, M.; Marek, M. Computational Methods in Bifurcation Theory and Dissipative Structures; Spring-Verlag: New York, 1983.
Received for review July 1, 2002 Revised manuscript received November 14, 2002 Accepted November 21, 2002 IE020493B
