HUNGARIAN JOURNAL OF INDUSTRIAL CHEMISTRY VESZPRÉM Vol. 33(1-2). pp. 49-56. (2005) STABILITY AND BIFURCATION OF CONTINUOUS COOLING CRYSTALLIZERS B.G. LAKATOS1and N. MOLDOVÁNYI2 1Department of Process Engineering, University of Veszprém, H-8201 Veszprém, P.O. Box 158, HUNGARY 2Honeywell Process Solutions, H-1139 Budapest Petneházy u. 2-4, HUNGARY Stability properties and dynamic behaviour of continuous cooling crystallizers are analysed using a detailed moment equation model. A causal loop diagram between the variables reveals that the roots of instabilities lay in the interactions of the autoinhibition generated negative feedback, a positive feedback between the four leading moments of crystal size and a varying polarity feedback between the temperature and the moments. The stability of steady states is analysed by eigenanalysis of the Jacobian matrix and using the Mikhailov criterion. Stability maps and bifurcation diagrams are presented in the planes of different pairs of system parameters. Keywords: Cooling crystallizer, moment equation model, causal loop diagram, stability, bifurcation Introduction Non-isothermal continuous crystallizers, extensively used in the chemical industry, usually are very sensitive to both the external and internal (parameter) disturb-ances. This is because of the highly nonlinear kinetics, many temperature-dependent parameters, which are, in turn, also nonlinear, and different feedbacks between the variables and elementary processes taking place in crystallization processes. All these properties, as well as the interactions between the kinetics, fluid dynamics and crystal size distribution may give rise to different complexities in both the steady state and dynamic beha-viour of continuous crystallizers, a deeper understand-ing of which is important in relation to both the crystal- lization process itself, and to the operation, control, and design of industrial crystallizers. Since the observations by Miller and Seaman [1] on composition and crystal size distribution oscillations in industrial crystallizers, a number of works have dealt with their stability and dynamic behaviour under iso-thermal conditions [1-16], but less attention has been paid to the problem of taking into account the thermal effects. Tavare et al. [17] studied the temperature mul-tiplicity and stability of cooling MSMPR crystallizers using linear temperature- dependence of the solubility. Melikhov et al. [18] shown that in a circulation vacuum crystallizer the dissolving zone may exhibit two steady states. Both works presented also stability criteria but used simplified models, and no primary nucleation and size-dependent growth were taken into account. The present study addresses the stability and bifurca- tion phenomena of continuous cooling MSMPR crystal- lizers by means of a detailed moment equation model. The stability is examined using the Mikhailov criterion, and the effects of primary and magma-dependent sec- ondary nucleation are analysed. Bifurcation diagrams and simulation results concerning the dynamic behaviour of crystallizers, as well as the causal loop diagram revealing the important negative and positive feedback loops of the system, responsible for the complex behaviour of crystal-lizers are presented. 50 Mathematical model Population balance model Consider a continuous crystallizer the schematic repre-sentation of which is presented in Fig.1. Let us assume that the following conditions are satis- fied: - the crystallizer may be seeded; Fig.1. Schematic representation of a continuous cooling crystallizer - the crystals can be characterized by a linear dimen-sion L; - all new crystals are formed at a nominal size so that we assume ; 0≥nL 0≈nL - crystal breakage and agglomeration are negligible; - no growth rate fluctuations occur; - the overall linear growth rate of crystals G is size- dependent and has the form of power law expression: (1) ( ) ( aLcckG gsg +−= 1 ) - the primary nucleation rate Bp is described by the Volmer model: ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −= s e pp c c k kB 2ln expε (2) where ε is the voidage of suspension expressed as (3) 31 µε vk−= and µ3 is the third of the ordinary moments of the population density function ,which are defined as (4) ∫ ∞ == 0 ,...3,2,1,0,),()( mdLtLnLt mmµ - the secondary nucleation rate Bb is described by the power law relation ( ) jbsbb cckB 3µ−= (5) where the coefficients kg, kp and kb are functions of temperature expressed as gpg RT E kk ,,,exp0 =⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −= ιιιι . (6) Under such assumptions the population balance mo- del of the crystallizer consists of the following balance equations. Vapour Volume balance of the crystal suspension qq dt dV in −= (7) Feed subject to the initial condition (7a) 0)0( VV = Cooling water out Population balance equation governing the crystal size dynamics 0,0, )()( >∞<<−=+ tLqnnq L Gn V t Vn inin∂ ∂ ∂ ∂ (8) subject to the initial and boundary conditions: (8a) nLLLnLn ≥= ),()0,( 0Product 0,),(),,(lim ≥+= → tBeBetLnLccG bbpps LL n (8b) Cooling water in 0,0),(lim ≥= ∞→ ttLn L (8c) where n(L,t)dL expresses the number of crystals having size in the range L to dL at time t, while ep and eb are bin- ary existence variables of the nucleation rates, by means of which the alternative variations of nucleation can be controlled. Naturally we have the constraint 1≥+ bp ee (9) Mass balance of solvent ( ) svsvininin sv cqcq dt cVd εε ε −= (10) with the initial condition (10a) 0)0( svsv cc = Mass balance of solute ])1([ ])1([ ])1([ c cinininin c cq cq dt VcVd ρεε ρεε ρεε −+− −−+= −+ (11) with the initial condition (11a) 0)0( cc = Energy balance for the crystal suspension takes form [ ] mccccsvsvc inccinsvinsvincinin h cccsvsv VRHTCcCcCq TCcCq dt ()(ε cC TTUaV TCVcCcCVd )(])1()([ ])1[ )( )1()( ∆ρεε ρε ρεε −+−++− −− +−−= −++ +++ (12) 51 with the initial condition (12a) 0)0( TT = where Rmc denotes the global rate of production of crystal mass in a unit volume of suspension: ( ) ( )[ ]32)1( 3 3 )(])1[( µµρρ ρµρε ε acckkR dt kd dt d R g sgVcc batch cv batch c mc +−== == − = − (13) Energy balance for the cooling medium ( ) ( ) )( hhhhinhinhhhhh TTUaVTqTqCdt TCVd −+−= ρ ρ (14) subject to the initial condition . (14a) 0)0( hh TT = The dependence of saturation concentration on the temperature is described by the expressions: ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −= T a aTcs 1 0 exp)( (15a) or . (15b) 2210)( TaTaaTcs ++= Often, it is useful to complete the set of balance equations (7)-(14) by the equation for the equilibrium saturation concentration given as dt dT dT dc dt dc ss = (16) subject to the initial condition . (16a) 00 )( ss cTc = Therefore, the state of the continuous cooling crys-tallizer at time t≥0 is given by the sextuple [V(t),c(t), csv(t),T(t),Th(t),n(.,t)], and its dynamics is described by the population balance model formed by the mixed set of partial and ordinary differential equations (7)-(14). The time evolution of this system occurs in the state space R5×N that is the Descartes product of the vector space R5 of volume, concentrations and temperatures, and the function space N of population density func-tions. Consideration of dynamical problems of crystal-lizers in this product space, however, seems to be quite complex so that hitherto only a few works have studied dynamic problems in this space directly. Controllability analysis was carried out by Semino and Ray [21], while Lakatos and Sapundzhiev [22] studied the global sta-bility of crystallizers via Lyapunov's Direct Method. In the present study, we concentrate on a reduced case, applying a finite- dimensional state space model based on the four leading moments of the crystal size instead of the distributed parameter system of Eqs (7)-(14). Moment equation model Since the overall crystal growth rate (1) is a linear func- tion of size L, the population balance equation (8) can be converted into an infinite set of recursive ordinary differential equations for the moments (4) of the pop- ulation density function which take the form ( ) bbppinin BeBeV q dt d ++−= 0,0 0 µµ µ (17) ( ) )()( 1, mmgsgminminm accmkV q dt d µµµµ µ +−+−= − m=0,1,2,3... (18) subject to the initial conditions .0)0( mm µµ = Since the set of Eqs (7)-(16) can be closed by means of the first four leading moments, we reduce the infinite set of Eqs (17)-(18) to the set of equations governing these moments. In this way, expressing the first derivat- ives of all the state variables, taking into account the tem-perature dependence of the parameters, and applying the scale factors g ctgV sskks 333 00 6: −−= , , , s , g ctgV sskks 222 01 6: −−= g ctgV sskks −−= 102 3: Vk=:3 }max{ 1 : in c c s = , { }in T T s max 1 := , t V q s s s =: where st and sV can be chosen arbitrary, we introduce the following set of dimensionless variables csycsymsxVsv csvcsvmmmV ===== ,,3,2,1,0,, µ qsftsTszTszcsy qthThTscs ===== ,,,, ξ and scaled parameters cc sρα = g c − Vap kD 6= abD , , , , assk tg −= 10β g ctgp sskk 343 00 −− bg ctgb j V sskkk −−−−= 34300 16 R sE Tg=gβ , R sE Tp=pβ , R sE Tb b =β , ht c Cs Uas =κ , h cT C Hs ∆ γ −= , , 00 asb c= T c s as b 11 = or , 11 asb T= 2 2 2 T c s as b = . Then, the dimensionless scaled equations in the state space form are given as: ff d dv in −=ξ (19) ( ) bbppinin eexxv f d dx ΘΘ ξ ++−= 00 0 (20) ( ) ( ) ( )101,11 exp xxzyyxxv f d dx gg sin in β β ξ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −−+−= ( ) ( ) ( )212,22 2exp xxzyyxxv f d dx gg sin in β β ξ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −−+−= (22) ( ) ( ) ( )323,33 3exp xxzyyxxv f d dx gg sin in β β ξ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −−+−= 52 (23) ( ) ( ) ( ) ( ) ( 32 3 3 ,3 3exp )1( )( 1 1 xxyy zx y yy xv xf d dy g s g in inin β βα ξ +−⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − − − − −− − − = ) (24) ( ) ( ) ( ) ( ) ( 32 3 , 3 3 3exp )1( )1 1 xxyy zx y yy xv xf d dy g s gsv svinsv ininsv β β ξ +−⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − − + +− − − = ) (25) ( ) ( ) ( 32 3exp )( xxyy z zzzz v f d dz g s g hin inin β β φ αγ φ κ φ φ ξ +−⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −+ +−−−= ) (26) ( ) )(, h hh hinh h hh zz v v zz v f d dz −+−= α κ ξ (27) ( ) ( ) ( ⎥ ⎥ ⎦ ⎤ +−⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −+ ⎢ ⎣ ⎡ −−−= 32 3exp )( xxyy z zzzz v f dz dy d dy g s g hin ininss β β φ αγ φ κ φ φ ξ ) (28) ( ) 331 xC C y C C y C C x h c h c sv h sv αφ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ +−= and ( ) in h c in h c insv h sv inin xC C y C C y C C x ,3,,31 αφ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ +−= (29) subject to the initial conditions 00 0000 )0(,)0( )0(,)0(,)0(,)0( hh svsvmm zzzz yyyyxxvv == ==== where ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −−= z y y k xD p s e app β Θ exp ln exp)1( 2 3 (31 a) and ( ) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −−= z xyyD bjbsabb β Θ exp3 . (31b) 3. 3) The causal loop diagram in Fig.2 exhibits an impor- tant positive feedback between the four leading moments x The moment equation model (19)-(29) of the con- tinuous cooling crystallizer is, in principle, a finite di- mensional dynamical system, generated by reducing the infinite dimensional population balance model without any simplifying assumptions. The reduced state is formed by the ninetuple ( 93210 ,,,,,,,, R∈= hsv zzyyxxxxvx ) of scaled, dimensionless variables and the time evolu- tion of this system occurs in a bounded region of the state space, defined in the following manner 9RΩ ⊂ ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ <≤< <≤≤≤≤ ≤≤≤≤ ≤≤≤≤≤≤ = max,min,max max33 67.0 max022 67.0 max011max00 max ,0 ,0.10,0 ,0,0 ,0,0,0 hhh svcsvin zzzzz xxxAx xAxxx syyyvv ρ xΩ thus the feasible region of solutions of the dynamical system (19)-(29) is formed by the compact domain (32) of the state space R9. Qualitative feedback analysis Since any instability problem can be viewed as a feed- back of properly phased signals, this difference can be explained qualitatively analysing the causal loop dia- gram, used extensively in system dynamics [23,24]. Since a positive feedback destabilizes, while negative feedbacks stabilize the systems, in principle, their inter- actions turn to be decisive for the dynamics of crystal- lizer. Because of the nonlinearities of crystallizer, the dominance of the negative and positive feedbacks may change in time during the course of the process what can generate diverse patterns of behaviour of the system. Let us consider the causal loop diagram of the mo- ment equation model in the basic size-independent case, presented in Fig.2, similarly to that derived by Lakatos [14] for an isothermal MSMPR crystallizer. Here, we show only the main negative and positive causal links between the variables of Eqs (19)-(28). From the point of view of dynamics, the following feedback loops are the most interesting parts of this causal loop diagram: 1. 1) There is a negative feedback in each of Eqs (19)- (28) due to the draw off of crystals and solution. These effects arise in all continuous flow systems naturally, and under certain conditions those by themselves are capable of stabilizing the crystallizer. 2. 2) There exists an important negative feedback loop y→x2→y generated by autoinhibition of the supersatura- tion: nucleation and growth of crystals relieves supersa- turation and inhibits further crystallization process. 0→x1→x2→x3→x0 that is closed through the nucleation rate. 4. 4) Depending on if the crystallization process is exo- thermic or endothermic there exists also a positive or ne- gative feedback between the temperature and the second moment of crystal size z→x2→z. In any case, however, increasing temperature affects the supersaturation neg- atively. Note that the concentration of solvent depends on a number of other variables but there is no feedback from ysv into the remaining variables, i.e. the concentration of solvent has no essential influence on the dynamic be- haviour of crystallizer while the volume of the suspen- sion has no effect on the dynamics at all. 53 The dynamic behaviour of crystallizers differs sig-nificantly on if the primary or secondary nucleation is the dominant mechanism of producing the new crystals [10, 11]. The main reason of this difference lays in the di-verse nature of positive feedbacks caused by the different forms of the nucleation rates. Namely, derivating both sides of Eq.(20) with respect to x3 we obtain inequality 00 3 <⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ξd dx dx d (33) for primary nucleation, i.e. in this case the positive feed- back is rate-decreasing. Similarly, for secondary nuclea- tion the inequality 00 3 >⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ξd dx dx d (34) holds, i.e. in this case the positive feedback is rate- increasing, forming a feedback of autocatalytic nature. Fig.2. Causal loop diagram of the basic size-independent case showing only the main feedback loops Stability and dynamic behaviour The qualitative analysis has shown that the variables v and ysv does not influence the dynamic behaviour sig-nificantly, while zh, being, in principle, a control var-iable affects the system only through modulating the temperature z. The behaviour of the crystallizer in the neighbour-hood of a steady state may be deduced by eigenanalysis of the Jacobian matrix of Eqs (19)-(28) at this state. The structure of the Jacobian matrix is shown in Fig.3 that supports the conclusions of the quality analysis. Determining the eigenvalues of the Jacobian matrix the characteristic polynomial of the system can be built up by means of which the Mikhailov plot is constructed and the Mikhailov stability criterion [25] can be applied. This criterion seems to be useful for analysing the stability of crystallizers, since it provides a simple gra-phical means for testing the stability of polynomials, al-lows a deep understanding and visual presentation of the phenomenon and can be formulated also in algebraic terms. ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = xxxxxxxx xx xxxxxxxx xxxxxxx xxxxxx xxxxxx xxxxxx xxxxxx xxxxxx 00 00000000 00 000 0000 0000 0000 0000 0000 0000000000 J Fig.3. The structure of the Jacobian matrix of system (19)-(28) The Mikhailov plots of the system (19)-(28) are shown in Fig.4 as a function of parameter ke in the case of primary nucleation, i.e. for eb=0. Since the Mikhailov stability criterion says that an nth order real polynomial is asymptotically stable if and only if the Mikhailov plot, starting at ω=0 from the positive real axis of the complex plane goes through n quadrants in turn, the plots in Fig.4 for ke=1.0 and ke=3.0 indicate stable states, while for ke=5.0 and ke=10.0 we have unstable states. The basic Negative feedback: autoinhibition Positive feedback: interaction of the moments Varying polarity feedback: exotermic-endotermic crystallization Negative feedback: withdrawal x1 x2 x3 y - - - - - x0 - + + + + + z - if γ>0, then + if γ<0, then - + + if γ>0, then + if γ<0, then - - 54 values of the parameters in all simulation runs were: Dap=5000, Dab=2000, v=1, fin=f=1, g=1, b=3, j=0.4, α=35, β=0, βg=1.1, βb=2.6, βp=3.45, b0=6, b1=3, b2=0.0,, φin=8.8, αh=30, vh=1, fh=1. Fig.4. Mikhailov plots of crystallizer for primary nucleation as a function of nucleation parameter ke These curves are shown only for smaller values of ω, since for large ω values the polynomial constructed on the basis of eigenvalues of the Jacobian matrix (19)-(28) behaves like (jω)9. Therefore, the behaviour at small ω values entirely determines the stability of the system. The Mikhailov plot can also be used to count the number of unstable zeroes, i.e. if it goes through only (n-k) quadrants in turn, then the polynomial has k un- stable zeros. Plots ke=5.0 and ke=10.0 go through quad-rants 1-2-3-2-3-4-1 in turn, what means that it goes through only 7 quadrants in turn thus these plots repre-sent unstable states showing two eigenvalues with pos-itive real parts. As a consequence, these eigenvalues are conjugate pairs so that these instabilities mean limit cycle oscillations. Indeed, Fig.5 shows the projection of development of the limit cycle oscillations for ke=10.0 into the R3 subspace of variables (x0,y,z). 1 2 3 0.750.76 0.770.78 0.790.8 0.81 0.94 0.96 0.98 1 1.02 Fig.5. Development of the limit cycle oscillations for ke=10.0 in the R 3 subspace of variables (x0,y,z). Fig.6 presents the bifurcation diagram ke–yS in which two Hopf bifurcations are seen at ke≈3.3 and ke≈15.9. Between these points the crystallizer generates limit cycle oscillations with varying amplitudes while outside this interval the crystallizer exhibits stable steady states. Fig.6. Bifurcation diagram ke-yS for primary nucleation with amplitudes of limit-cycle oscillations When the secondary nucleation becomes the dom- inant mechanism of forming the new crystals, i.e. when ep=0, then the crystallizer exhibits steady state multipli- city and may have one, two, or even three steady states depending on the values of parameters. The detailed multiplicity analysis of nonisothermal crystallizers will be presented elsewhere, but it can be proved easily that for j>0 the system (19)-(28) always has a trivial boundary (washout) steady state. Under such conditions the crystallizer can not be ignited because of the small nucleation rate so that the washout steady states are always stable. Fig.7 presents the bifurcation diagram j-yS in the case of magma-dependent secondary nucleation where two regions are set by a boundary bifurcation point at j≈0.8. In the first region the crystallizer has two stable steady states, while in the second region only a unique washout steady state exists. Fig.7. Bifurcation diagram j-yS for secondary nucleation with boundary bifurcation of the washout steady state The bifurcation diagram Dab-yS of similar form is shown in Fig.8 in the case of magma-dependent second- ary nucleation where two regions are set by a boundary bifurcation point at Dab≈460. Here, the first region ex- hibits a unique washout steady state while two stable steady states exist in the second region. yS Im Re ke=1.0 ke=3.0 ke=5.0 ke=10.0 ω=0 ke=3 ke=15 ke Boundary bifurcation x0 y z ke=10 Two steady states Unique washout steady state yS j 55 Fig.8. Bifurcation diagram Dab-yS for secondary nucleation with boundary bifurcation of the washout steady state In the case of magma-dependent secondary nuclea-tion limit cycle oscillations have not been observed in the feasible region of parameters. Fig.9 presents some tran-sients of dimensionless concentration as a function of parameter Dab. Increasing Dab increases the willingness of crystallizer to generate damped oscillations in the tran- sient behaviour but next the crystallizer achieves stable stationary states. This phenomenon is important when designing the controlling system for s crystallizer. Fig.9. Damped oscillations of transients of crystallizer in the case of magma-dependent secondary nucleation Conclusions The stability and bifurcation analysis of continuous cooling crystallizers, carried out by means of a detailed moment equation model, has shown that there is a signi-ficant difference between the dynamic behaviour of crystallizers with primary and secondary nucleation. Continuous cooling crystallizers may exhibit limit cycle oscillations and multiple steady states. In the case of magma- dependent nucleation there exists always a washout steady state where the crystallizer is not ig-nited. The Mikhailov plots and criterion, constructed on the basis of eigenanalysis of the Jacobian matrix al-lowed a very useful visual form of stability analysis in the frequency domain. The causal loop diagram showing the negative and positive causal links between the model variables of crystallizer revealed that a positive feedback loop exists between the four leading moments of crystal size dis- tribution closed through the nucleation rate. The auto- inhibition generated negative feedback is closed mainly through the crystal growth rate although during the onset of crystallization the nucleation rate proves to be the dominant factor also in this feedback loop. The tem- perature with the second moment and solute concentra- tion forms feedbacks of positive or negative polarity de- pending on if the crystallization process is exothermic or endothermic. Interactions of these feedbacks bet-ween the variables with the kinetic nonlinearities appear to be decisive for the instabilities of crystallizers. Simulation studies concerned with the steady state multiplicity and stability patterns, as well as with the dynamic behaviour of cooling crystallizers allow us to conclude that continuous crystallizers have a broad range of complex behaviour in both the steady and dy-namic states. A further examination of these phenomena seems to be useful for both a better understanding of the crystallization process itself, and for improving the ope- ration, control, and design methods of industrial crys- tallizers. Acknowledgement This work was supported by the Hungarian Research Foundation under Grant T034406. Notation a - constant of the crystal growth rate [m-1] b - exponent of secondary nucleation rate Bp - primary nucleation rate [no m -3s-1] Bb - secondary nucleation rate [No m -3s-1] c - concentration of solute [kgm-3] cs - equilibrium saturation concentration [kgm -3] Dap - dimensionless parameter for primary nucleation Dab - dimensionless parameter for secondary nucleation g - exponent of crystal growth rate G - crystal growth rate [ms-1] Im - imaginary part of a complex number j - exponent of secondary nucleation rate ke - parameter of primary nucleation rate kg - rate coefficient of crystal growth [m 3g+1 kg-g s-1] kp - rate coefficient of primary nucleation [no m -3s-1] kb - rate coefficient of secondary nucleation [no m 3b-3 kg-b s-1] kV - volume shape factor L - linear size of crystals [m] n - population density function [no m-4] Re - real part of a complex number sc - scale factor of the concentration [kg -1m3] Boundary bifurcation yS Two steady states Unique washout steady state Dab y Dab 500 1000 5000 10000 ξ 56 sm - scale factor of the m th order moment of n (m=0,1,2,...) xm - m th order dimensionless moment (m=0,1,2,...) y - dimensionless concentration of solute Greek letters α - dimensionless parameter ß - dimensionless parameter γ - dimensionles crystal growth rate ε - viodage of suspension µm - m th order moment of n [mm-3] Θ - dimensionless nucleation rate ρc - density of crystals [kgm -3] ξ - dimensionless time Subscripts 0 - initial value in - inlet value p - primary nucleation b - secondary nucleation S - steady state References 1. MILLER, P. and SAEMAN, W.C.: Chemical Eng- neering Progress, 1947, 43, 667-675. 2. RANDOLPH, A.D. and LARSON, M.A.: AIChE Journal, 1962, 8, 639-645 3. NYVLT, J. and MULLIN, J.W.: Chemical Engineer- ing Science, 1970, 25, 131-147. 4. YU, K.M. and DOUGLAS, J.M.: AIChE Journal, 1975, 21, 917-923. 5. SONG, Y.H. and DOUGLAS, J.M.: AIChE Journal, 1975, 21, 924-925. 6. JERAULD, G.R., VASATIS, Y. and DOHERTY, M.F.: Chemical Engineering Science. 1983, 38, 1675- 1681. 7. WITKOWSKI, W.R. and RAWLINGS, J.B.: Proceed- ings of the American Control Conference, Min- neapolis, MN, 1987, 1400. 8. RANDOLPH, A.D. and LARSON, M.A.: Theory of particulate processes. Second Edition. Academic Press, New York, 1988. 9. BUEVICH, YU.A., MANSUROV, V.V. and NATA- LUKHA, I.A.: Chemical Engineering Science, 1991, 46, 2573-2578. 10. LAKATOS, B.G.: Computers and Chemical Engin- eering, 1994, 18, S427-S431. 11. LAKATOS, B.G. and SAPUNDZHIEV, TS.J.: ACH – Models in Chemistry, 1995, 132, 379-394. 12. KIND, M. and NIEKEN, U.: Chemical Engineering and Processing, 1995, 34, 323-328. 13. PATHATH, P.K. and KIENLE, A.: Chemical Engin- eering Science, 2002, 57, 4391-4399. 14. LAKATOS, B.G.: Chemical Engineering Transac- tions, 2002, 1 (3), 1149-1154. 15. MOTZ, S., MITROVIC, A., GILLES, E.D., VOLLMER, U. and RAISCH, J.: Chemical Engineering Science, 2003, 58, 3473-3488. 16. YIN, Q.X., SONG, Y.M. and WANG, J.K.: Industrial and Engineering Chemistry Research, 2003, 42, 630-635. 17. TAVARE, N.S., GARSIDE, J. and AKOGLU, A.: AIChE Journal, 1985, 31, 1128-1135. 18. MELIKHOV, I.V., ZELENKO, V.L., PODKOPOV, V. M. and IVAKIN, I.V.: Teoreticheskie Osnovy Khi- micheskoj Tekhnologii, 1992, 26, 56-63. 19. SEMINO, D. and RAY, H.: Chemical Engineering Science, 1995, 50, 1805-1824. 20. LAKATOS, B.G. and SAPUNDZHIEV, TS.J.: Bulga- rian Chemical Communications, 1997, 29, 28-38. 21. KIRKWOOD, C.W.: System Dynamics Methods. A Quick Introduction. Arizona State University, Tempe, 1988. 22. RICHARDSON, G.P.: System Dynamics Review, 1995, 11, 67-88. 23. MACFARLANE, A.G.J.: Frequency-Response Me- thods in Control Systems. IEEE Press, New York, 1999.