Microsoft Word - 01_R.doc HUNGARIAN JOURNAL OF INDUSTRIAL CHEMISTRY VESZPRÉM Vol. 35. pp. 7-17 (2007) POPULATION BALANCE MODELLING OF CRYSTALLISATION PROCESSES B. G. LAKATOS University of Pannonia, Institute of Chemical and Process Engineering, Department of Process Engineering, H-8201 Veszprém, P.O. Box 158, HUNGARY The work presents and discusses the basics of the population balance approach and multi-dimensional population balance models for describing different aspects of solution crystallisation. A two dimensional model is presented for describing evolution of needle-shape crystals characterised by two size dimensions in cooling crystallisation, and a two-population model is developed for describing micromixing of solution, and its effects on the size distribution of crystals and their aggregates in reaction crystallization. Moment equation reductions for joint moments of internal variables are developed in both cases which are closed by means of cumulant neglect closure models. The properties and behaviour of crystallization processes described by the models are investigated by simulation. Keywords: Cooling and reaction crystallisation, multivariable population balance, needle-form crystals, micromixing Introduction In Model-Driven Chemical Engineering which refers to the systematic use of models as primary engineering artefacts throughout the engineering lifecycle models are considered as first class entities, and due to their cognitive and methodological values serve as fundamental tools of understanding, analysis, design, operation, control and managing of real engineering processes and processing systems. This technique states that any specification should be expressed by models, which are both human and machine understandable. The Department of Process Engineering attempts to organize new efforts in these directions by proposing a framework (1) to clearly define methodologies, (2) to develop models at any level of abstraction, and (3) to organize and automate the testing and validation activities. Models, depending on what they represent, can reside at any level of abstraction, and can be restricted to address only certain aspects of the system. Well defined, adequate mathematical models usually are of great predictive power and together with the corresponding computer models have become efficient and useful aspects of engineering, especially treating complex engineering processes and processing systems. In chemical and process engineering such complex systems are, among others, those multiphase processes in which at least one of the phases take part in dispersed state, i.e. in form of solid particles, liquid drops or gas bubbles. These systems can not be modelled adequately by means of the classical state variables only since the characteristic properties of the dispersed elements, among others their size, size distribution, habit, shape, structure, or configuration in continuous phase, have to be taken into consideration as well. For these reasons, population balance models in which population balance equations describe the properties and behaviour of dispersed elements, appear to be adequate mathematical models of disperse systems. In chemical engineering, Hulburt and Katz [1] formulated the population balance equation using the ideas of statistical mechanics, while Randolph and Larson [2] based their formulation on continuous mechanics aiming to describe crystallization processes. These formulations have found a lot of applications in modelling particulate systems [2-5], systems containing liquid drops [6-8], as well as gas bubbles [9]. Ramkrishna [10] presented an excellent state of the art of this research field. While taking part in technological processes disperse elements interact intensively with the continuous phases, and also with the counterparts of the dispersed one. As regards the dispersed phase, here two types of interactions can be distinguished: 1) interactions causing changes of the dispersed elements themselves, and 2) interactions inducing changes of entities carried by the dispersed elements. Continuous growth and shrinkage, breakage and aggregation-agglomeration or coalescence of dispersed elements can be mentioned in the first group, while exchange of momentum, mass exchange of chemical species, exchange of heat and electrical charge may be included into the second group. In effective modelling of disperse systems, these interactions and processes induced by those should to be taken into consideration. When modelling solid dispersed phase, i.e. particulate systems, in most part only those inter-particle interactions have been accounted for which cause changes of the particles themselves, i.e. their breakage, aggregation 8 and agglomeration. Recently, developing a new, Markov process-based formulation of the population balance equation [11, 12], direct inter-particle heat exchange processes [13-16] and, mostly for modelling direct interactions between fluid particles, also mass exchange processes [17-20] have been included into the population balance models, extending in this way the range of effect of the population balance modelling approach. The aim of the present work is to present and discuss the basics of the population balance approach and multi- dimensional population balance models for describing different aspects of crystallisation from solution. A two dimensional model is developed for describing evolution of needle-form crystals characterised by two size dimensions in cooling crystallisation, and a two- population model is presented for describing the mixing effects of solution on the size distribution of crystals and their aggregates in reaction crystallization. Moment equation reductions for joint moments of internal variables are developed in both cases which, when the sets of equations are unclosed, are closed by means of cumulant neglect closure models. The properties and behaviour of the systems described by the models are investigated by simulation. Two-dimensional population balance models Population balance model of a cooling crystallizer for needle-shaped crystals For describing crystallization from solution, which is a widely used cleaning, separation and particle producing technique in the chemical, petrochemical and pharmaceutical industry, population balance models have been proved adequate modelling tools. In this enormously important unit operation, the crystal population is the main product thus the goal of operating crystallizers is to produce crystals of prescribed form, habit, size, and size distribution. The task of describing this process is to predict the properties and behaviour of the whole crystal population from the behaviour of single crystals interacting with the solution and crystals environment. The population is modelled by the density of suitable chosen variables as models of crystals’ characteristic properties variation of which is described by the population balance equation. Most of the modelling studies of crystallisation use one-dimensional characterisation of crystals. However, crystalline particles often, especially in the pharmaceutical industry, exhibit much more complex habit, as it is illustrated in Fig. 1. Crystals in Fig. 1b and 1c can be described adequately only by multi-dimensional models, characterising the crystal and crystal growth separately for the main crystallographic faces. Here we consider needle-shaped crystals, produced in a cooling crystallizer shown schematically in Fig. 2. Such needle-shaped crystals can be characterised by two size dimensions L1 and L2, while the crystal population is described by the population density function n(., ., t) by means of which n(L1, L2, t)dL1dL2 expresses the number of crystals from the size domain (L1, L1 + dL1) × (L2, L2 + dL2) in a unit volume of suspension at time t. Since the volume of crystals is an important property a volumetric shape factor kV is also defined by means of which the crystal volume can be expressed using only these two size dimensions even in case of different shapes as it is shown in Fig. 1c. Figure 1: Three different crystal systems Figure 2: Continuous cooling crystallizer Let us now assume that: - the working volume of the crystallizer is constant; - all new crystals are formed at a nominal size L1, n ≅ L2, n ≅ Ln ≥ 0, so that we can assume: Ln ≈ 0; - crystal breakage and agglomeration are negligible; - the overall linear growth rates of the two habit faces G1 and G2 are size independent, and have the form of power law expressions: G1 = kg1(c – cs) g1 (1) and G2 = kg2(c – cs) g2 (2) - the primary nucleation rate Bp is described by Volmer’s model Feed Product Vapour Cooling water out Cooling water in q, cin, csv,in,Tin q, c, csv T, n qh, Th,in qh, Th L1 L3 L2 L1=L2=L3 L3 L1 L2 L1≠L2>>L3 L3 L1 L2 L1=L2=∞<< −=++ tiL nn V q L nG L nG t n i in∂ ∂ ∂ ∂ ∂ ∂ (9) subject to the initial and boundary conditions: n(L, 0) = n0(L), L ≥ Ln (10a) 0,),(),,(lim ≥+= → tBeBetLnLccG bbpps LL n (10b) 0,0),(lim ≥= ∞→ ttLn L (10c) where ep and eb are binary existence variables of the nucleation rates, by means of which the alternative variations of nucleation can be controlled. Naturally we have the constraint ep + eb ≥ 1 (11) Mass balance of solvent: ( ) ( )svsvininsv ccV q dt cd εε ε −= (12) with the initial condition csv(0) = csv0. The mass balance equation of solute: ])1([ ])1([ ])1([ c cininin c c V q c V q dt cd ρεε ρεε ρεε −+− −−+= −+ (13) subject to the initial condition c(0) = c0. The energy balance equation for the crystal suspension takes form [ ] TCcCcC RH q TCcCcC V q RHTTUa dt TCcCcCd ccsvsvc mcc inccinsvinsvincin mcch cccsvsv ])1()([ )( ])1()([ )()( )1()( ρεε ρεε ρεε −++ − − −+++ +−+−−= = −++ Δ Δ (14) subject to the initial condition T(0) = T0. Here, Rmc denotes the global rate of production of crystal mass in a unit volume of suspension. Energy balance for the cooling medium: ( ) ( ) )( hhhin h h hh hh TTUaVTT V q C dt TCd −+−= ρ ρ (15) subject to the initial condition Th(0) = Th0. The dependence of saturation concentration on the temperature is described by the expression: cs(T) = a0 + a1T + a2T 2. (16) Often, it is useful to complete the set of balance equations (9)-(15) by the equation for the equilibrium saturation concentration given as ( ) dt dT Taa dt dT dT dc dt dc ss 21 2+== (17) subject to the initial condition cs(T0) = a0 + a1T0 + a2T0 2. Therefore, the state of the crystallizer is given by the quintuple [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 (9)-(15). The time evolution of this system occurs in the state space R4×N that is the Descartes product of the vector space R4 of concentrations and temperatures, as well as the function space N of population density functions. In the present study, we concentrate on a reduced case, approximating the distributed parameter system of Eqs (9)-(15) by a finite- dimensional state space model using the six leading joint moments of crystal sizes. 10 Moment equation reduction and scaling The moment equation reduction of the population balance model (9)-(15) can be written for the joint moments of crystal sizes (5) by means of which the mean values of sizes are given as )( )( )( 0,0 0,1 1 t t tL μ μ = and )( )( )( 0,0 1,0 2 t t tL μ μ = (18) The infinite hierarchy of the moment equations corresponding to Eq. (9) takes the following form ( ) bbppin BeBeV q dt d ++−= 0,0,0,0 0,0 μμ μ (19) ( ) ...3,2,1,,1,2 ,11,,, , =+ ++−= − − mkmG kG V q dt d mk mkmkinmk mk μ μμμ μ (20) This infinite set of equations can be closed at any order but natural closing occurs at the third order joint moment μ1, 2 since, because of Eq. (6), this moment is required for the mass balance of solute (13). To do that equations for moments μ0,0 μ1,0 μ0,1 μ0,2 μ1,1 are also required. Since the system of moment equations often become stiff it is advisable to perform scaling that by defining useful scale parameters and dimensionless scaled variables to allow controlling conditioning of the system. Introducing the scale parameters s1, 2 = kV, s0, 2 = kVψ1, s1, 1 = 2kVψ2, s1, 0 = s1, 1ψ2 = 2kVψ2ψ2, s0, 1 = s1, 1ψ1 = 2kVψ2ψ2 s0, 0 = 2kV ψ1ψ2ψ2, sc = 1/cin, sT = 1/Tin and st, where 1 1 1 0,1 g ctg ssk −−=ψ , 2 2 1 0,2 g ctg ssk −−=ψ , the dimensionless scaled variables xk, m = sk, m μk, m, k, m = 0, 1, 2, 3, ysv = sc csv, y = sc c, ys = sc cs, z = sTT, zh = sT Th, ξ = st t and scaled parameters cc sρα =: , 2 21 1 02: ψψ −= tpVap skkD , 2 21 1 0 12: ψψbctb j Vab sskkD −−−= , gbp R sE T ,,,: == ιβ ιι ht c Cs Uas :=κ , )(: c h T H C s Δγ −= 00 : asb c= , T c s as b 11 := 2 2 2 : T c s as b = . the closed set of moment equations takes the form ( ) bbppin eexxd dx ΘΘ ξξ ++−= 0,0,0,0 0,0 1 (21) ( ) 0,00,1,0,10,1 1 1 xxx d dx gin Θξξ +−= (22) ( ) 0,01,0,1,01,0 2 1 xxx d dx gin Θξξ +−= (23) ( ) 0,11,01,1,1,11,1 21 1 xxxx d dx ggin ΘΘξξ ++−= (24) ( ) 1,02,0,2,02,0 2 1 xxx d dx gin Θξξ +−= (25) ( ) 2,12,1,2,12,1 1 Ψ ξξ +−= xx d dx in (26) ( ) ( ) ( ) 2,1 33 ,3 )1( )( 1 11 Ψ α ξξ x y yy x x d dy in in − − −− − − = (27) ( ) ( ) ( ) 2,1 3 , 3 3 )1()1 11 Ψ ξξ x y yy x x d dy sv svinsv insv − +− − − = (28) ( ) 2,1)( 1 Ψ φ αγ φ κ φ φ ξξ +−−−= hin in zzzz d dz (29) ( ) )(1 , h hh hinh h h zz v v zz d dz −+−= α κ ξξ (30) subject to the initial conditions xk, m(0) = xk, m0, k, m = 0,0; 0,1; 1,0; 1,1; 0,2; 1,2 ysv(0) = ysv0, y(0) = y0, z(0) = z0, zh(0) = zh0. In Eqs (21-30) ( ) 1,12,01,12,02,1 21, xxxx gg ΘΘΨ += (31) ( ) 2,12,11 xC C y C C y C C x h c h c sv h sv αφ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ +−= (32) ( ) in h c in h c insv h sv inin xC C y C C y C C x ,2,1,,2,11 αφ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ +−= (33) ( ) 2,1,exp =⎟⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ −−= i z yy ii i gg sg β Θ (34) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −−= z y y k xD p s e app β Θ exp ln exp)1( 2 3 (35) and ( ) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −−= z xyyD bjbsabb β Θ exp3 . (36) According to Eq. (17), the set of Eqs (21-30) can be completed using the equation describing the variation of the dimensionless saturation concentration 11 ( ) ( ) ( ) ( )1,12,02,121 21 ,2 )( 1 2 xxzbb zzzzzbb d dy hin ins Ψ φ αγ φ κ φ φ ξξ ++ +⎥ ⎦ ⎤ ⎢ ⎣ ⎡ −−−+= (37) Note that when L1 = L2 = L, we can write μk,m = μk+m, k, m = 0, 1, 2, 3, ..., and the set of equations (20)-(30) reduces to the one dimensional moment equation model of cooling crystallizers. Simulation results and discussion The system of ODEs (21)-(30) was solved in MATLAB environment using the process and kinetic parameters presented in Table 1. Kinetic parameters were chosen followed the works by Ma et al. [27] and Briesen [28]. Simulation runs have been carried out to investigate the influence of nucleation on the crystal shape under the same crystal growth and operation conditions. Table 1: Process and kinetic parameters used in simulation of cooling crystallizer V = 10 m3 kp0 = 1.6·10 30 g1 = 1.48 qin = 10 -3 m3 s-1 ke = 2.1 Eb = 9.0·10 4 cin = 350 kg m 3 Ep = 5.0·10 4 mol kg1 = 12.2·10 -6 m s-1 Tin =90 °C kb0 = 1.0·10 16 g2 = 1.75 Th = 20 °C b = 3.0 j = 1.5 kg2 = 10.08·10 -7 m s-1 a0 = 0.2087 a1 = –9.76·10 -5 a2 = –9.30·10 -5 Fig. 3 presents the effects of parameter ke of primary nucleation on the steady state mean size values 〈L1〉 and 〈L2〉. The crystallizer produced practically the same mass of crystals in each run, i.e. about 90 kg m-3 of suspension. Figure 3: Effects of parameter ke of primary nucleation on the steady state mean size of the crystallographic faces Significantly different diagrams were obtained when secondary nucleation was the dominant crystal producing process. As it is shown in Fig. 4, in this case crystals were able to achieve size about 1 mm in length while producing almost the same amount of crystalline mass as in the case of primary nucleation. Naturally, the differences between the numbers of crystals took some orders of magnitude. The operation mode used in simulation assured good cooling conditions as it is illustrated in Figs 5 and 6. Projections of trajectories of crystallizer into the 3-D Figure 4: Effects of parameter b of secondary nucleation on the steady state mean size of the crystallographic faces T–c–μ0 subspace are presented in Fig. 5 for primary nucleation and in Fig. 6 for secondary one. It is seen well that the transients differ significantly depending on the nature of nucleation. Figure 5: Trajectories of the crystallizer in the 3-D T-c- μ0, 0 subspace in the case of primary nucleation b 1L 1L 2L 2L ke 1L 2L 1L 2L c T 12 Figure 6: Trajectories of the crystallizer in the 3-D T-c- μ0,0 subspace in the case of secondary nucleation Two population models Generalized coalescence/redispersion model for micromixing The physical phenomena occurring in solution crystallizers are very complex, involving the interactions between large-scale fluid mixing, fine-scale micromixing and complex crystallization kinetics as well as, in reaction crystallizers, chemical reaction kinetics. The macroscopic models of the fluid elements adequately describe macromixing but a quantitative description of the influence of micromixing presents a more challenging problem. Most of the works dealing with the influence of micromixing on crystallisation processes was focused on reaction precipitation [21-24]. In crystallisation, however, fine-scale fluid mixing may also affect the crystal size distribution because of the mixing of fresh feed stream with the bulk solution of crystallizer. Then, direct interactions between the supersaturation inhomogeneities coupled with highly nonlinear nucleation and crystal growth kinetics may influence the process significantly, as it was demonstrated by means of a probability density function model [25], and, recently, by developing two- population models with the usage of coalescence- redispersion micromixing closures [17-19,26]. Here, an improved two-population model is presented in which the generalised coalescence-redispersion micromixing model [20] is utilized for describing micromixing of solution. Let us now assume the in a continuous crystallizer the supersaturation is generated by a chemical reaction according to the scheme ↑⎯→⎯+ CBA k . Further, we assume that the fluid phase in the crystallizer can be visualised as consisting of a large number of fluid ele ments in which concentrations and supersaturation can be treated as deterministic quantities. The fluid elements, moving stochastically in the crystallizer, interact intensively with the counterparts and the existing crystals. Fluid elements exhibiting different concentrations of species A, B and C change some amounts of mass between each other during their coalescence-redispersion interactions, inducing in this way equalization of concentration differences. The population of fluid elements is characterised by the density function p(.,.), where p(c, t)dc expresses the number of fluid elements having concentration in domain (c, c+dc) in a unit volume of suspension at time t. Here, c = (ca, cb, cc) denotes the concentration vector of species A, B and C. The density function p(.,.) describes, in essence, the distribution of concentrations of species in the fluid elements so that the average concentration observable on macrolevel is expressed as ∫ == m cbakdtpc t tc kk c cc 00 ,,,),( )( 1 )( π (38) where ∫= m dtpt c cc 0 0 ),()(π (39) provides the total number of fluid elements in a unit volume of suspension. In Eq. (38) expression c c d t tp t )( ),( )( 0π =Prob (40) is interpreted as the probability of there exists a fluid element having concentration in domain (c, c+dc) in a unit volume of suspension at time t. In the crystallizer, between the stochastically moving fluid elements and crystals the following interactions are assumed to occur: 1) fluid-fluid element, 2) fluid element-crystal, 3) crystal-crystal interactions. Interactions between the fluid elements are manifest as micromixing of solution that is described by the generalized coalescence-redispersion model [20]: [ ] ( ) ( ) ∫ ∫∫∏ ∫ ∫∫∏ = = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ −⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + − + + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ −⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + − −= X X X X iM 1 00 1 00 "'),"(),'()("' )'(212 "'),"(),'()(" )'(212 3 1 3 3 1 3 ccccω ccccω ω ω ddtptpdFcc cctS ddtptpdFcc cctS p k kk k kk c k col k kk k kk c k col k k ω δ ωπ ω δ ωπ (41) Here, ωk ∈ [0,1] is a random parameter with distribution function Fωk, characterizing the exchange process of the kth, k = a, b, c component of c between the fluid elements, while Scol stands for the frequency of binary collisions. When Fωk = 1(ωk), i.e. the mean value of ωk is zero, then no mixing of the k th species occurs on microlevel, but Fωk = 1(ωk – 1), i.e. the mean value of ωk equals 1, denotes maximal micromixing under the given macrolevel conditions. μ 0,0c T 13 As a consequence, beside the crystal population the solution in crystallizer is treated as a second population of the fluid elements, forming in this way a two- population system of interacting populations of two different dispersed elements. Crystallization kinetics of under non-perfect micromixing conditions In the crystallizer, simultaneously with the chemical reaction of species A and B producing the precipitating species C, nucleation of crystals and their subsequent growth occur, due to the direct fluid element-fluid element and fluid element-crystal interactions, whilst the possible crystal-crystal collision interactions may induce also aggregation of crystals. Under such conditions, it seems to be better to characterise crystals on the volume coordinate v, and the population of crystals by the population density function n(., t) by means of which n(v, t)dv expresses the number of crystals from the crystal volume interval (v, v+dv) in a unit volume of suspension at time t. As a consequence, now the motion of the two populations occur in a four-dimensional space of internal properties from which the first one is the volume coordinate of crystals, while the remaining three internal properties are formed by the concentration coordinates of fluid elements. Assuming that - all new crystals are formed with a nominal volume vn ≈ 0 and vn > 0, then the kinetic expressions for nucleation, growth and aggregation of crystals are as follows. Primary nucleation is considered as a result of fluid element interactions while secondary nucleation is a result of those between the fluid elements and crystals. If the intrinsic rate of nucleation is given by expression Bt(cc, cs, v) its macrolevel rate is expressed as ( ) ( ) bpdvdtvntp vccBtB nv sc ,,),(),( ,, 1 0 =× ×= ∫ ∫ ∞ ι μπ ιι cc mc 00 (42) In similar way, since crystals moving stochastically in the vessel are surrounded and suffer contacts with fluid elements of different composition, the macrolevel crystal growth rate is also expressed as ( ) ( )∫== mc cc 00 dtpvccG t tG dt dv sc ),(,,)( 1 π (43) Finally, the rate of aggregation of crystals by means of the frequency of collisions and the aggregation efficiency is given as ( ) ( ) ( ) ( )∫∫ ∫∫ −−+ +−−−= ∞ m m v asceff asceffa dtpdvtvntvvnvvvSccS dtpdvtvvntvnvvvSccStB c c cc cc 0 00 0 00 ),('),'(),'(,',', 2 1 ),('),'(),(,',, 1 )( ε πμ ε πμ 0 0 (44) Two-population model of reaction crystallization If the reaction crystallizer 1) is operated under isothermal conditions, 2) the frequency of coalescence-redispersion events is constant, and 3) the reaction between species A and B is quasi-linear the rate of which is given as Ra = –kcacb, Rc = kcacb, then the population balance equation for the fluid elements takes the form [ ] [ ] ( )[ ] ( ) ( ) ∫ ∫∫∏= ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ −⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + − +−= ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ −=− ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ X X 1 00 "'),"(),'()("' )'(212 ),( ),(),(),(),(),( ),( 3 1 3 ccccωc ccccc c ω ddtptpdFcc cc t tS tpS tptp V q tpRR c tpR c tpR ct tp k kk k kk c k col col in in c c b b a a k ω δ ωπ ε ε (45) In Eq. (45), expression 〈R〉 represents the rate of decrease of concentration of solute in a fluid element due to crystallization which, because of the rate expressions (42) and (43), can be expressed as an average over all crystals being in its environment: n ccc vB c G c R dt cd ιε ρ μ ε ρ − − − −== 0 (46) Boundary conditions for concentration variables specify closed system along the concentration coordinates except the input where any combination of feasible conconcentration values can appear in the inlet density function. This means that for any k ∈ {a, b, c} and any values of the remaining concentration variables 0),()( 0 = =kck tpR cc and 0),()( , = = mkk cck tpR cc . The second population balance equation of the model describes the crystals is formulated here as 14 ( ) ( ) ( )[ ] [ ] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )∫ ∫ ∫ ∫ −−+ +−− −−+−=+ ∞ m n m n v v sceffa v sceffa nin ddvtptvntvvnccSvvvS tt ddvtptvnccSvvvS tt tvn vvtBtvntvn V q v tvnvtG t tvn c 0 c 0 cc cc 00 00 '),(,',',',' 2 1 '),(,',,' , ),(),( ,,, πμ πμ δ ∂ ∂ ∂ ∂ (47) The boundary and initial conditions to Eq. (47) can be written, respectively, as follows: ( ) ( ) 0,,lim = +→ tvnvtG nvv and ( ) 0,lim = ∞→ tvn v (48) and p(c, 0) = p0(c) and n(v, 0) = n0(v) (49) Writing Eq. (47) it was assumed that the volume vn of nuclei is small but vn > 0, thus the nucleation rate here arises in Eq. (47) as a source term inside the support interval [0, ∞) of crystal volumes. The second consequence of this assumption is the first of the boundary conditions (48). Moment equation model The two-population model of the reaction crystallizer is formed by the population balance equations (45) and (47), respectively, for the fluid elements and crystals. Eq. (45), however, is a three-dimensional equation by itself so that a moment equation reduction requires the joint moments of concentrations, defined as ....2,1,0,,,),()(,, == ∫ nlmdtpccct m n c l b m anlm c cc 0 π (50) The infinite hierarchy of the moment equations corresponding to Eq. (45) takes the form ...2,1,0,,,)()( 1 )()( )()()()()()( )( 0 0 0 ,,,,,,, 0,0,0 ,,,, ,,,1,,1,1,11,,11,1,1 ,, =⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ −− −=+−++ ∂ ∂ ∑∑∑ = = = −−− −−++++ nlmttcbatSt V q t V q ttRntnktlktmk t t n k l j m i knjlimkjinkljminlmcolnlm innlm in nlmnlmnlmnlm nlm ππ π ππ π ε ε ππππ π (51) where the coefficients of micromixing terms are expressed as ∫ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −⎟ ⎠ ⎞ ⎜ ⎝ ⎛= 1 0 ,, )(22 1 aa im a i a mi dFi ma ω ωω ω (52a) ∫ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −⎟ ⎠ ⎞ ⎜ ⎝ ⎛= 1 0 ,, )(22 1 bb jl b j b lj dFj lb ω ωω ω (52b) ∫ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −⎟ ⎠ ⎞ ⎜ ⎝ ⎛= 1 0 ,, )(22 1 cc kn cc mi dFkk nc ω ωω ω (52c) Eqs (51) and (52) show that the micromixing terms are always closed so that the joint moment equations may be unclosed only due to the reaction terms as it seen in Eq. (51). As a consequence, the mean values of concentrations (38) are written as π1, 0, 0/π0, π0, 1, 0/π0 and π0, 0, 1/π0. The second part of the moment equation reduction of the two-population model of reaction crystallizer is based on the moments of the crystal volume, given as ( ) ( ) ...2,1,0, == ∫ ∞ kdvtvnvt nv k kμ (53) Assuming now hat aggregation of crystals is characterised by sum kernel in which the efficiency is equal 1, then Sa(v, v')Seff(cc, cs) = Sa·(v + v') (54) and the resulting hierarchy of moment equations takes the form ( ) ( ) [ ] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ...2,1,0,1 2 )()()( 1 0 1 0 0 101 ,1 =⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + + + − −+−=− ∑ += = −+ + − ktt j k t S t tttt S vtBtt V q tGk t t kj j jkj a kk a k nkinkk k μμ μ μ μμμμ μμμ ∂ ∂μ (55) Note that the zero order moment of the set (53) provides the total number of crystals, the first order moment expresses the total volume of crystals in a unit volume of suspension, while the second order one makes possible of predicting the dispersion of the crystal population. In this case the first three moments seem to be sufficient to characterise the crystal population. As a consequence, taking into account Eq. (51), the order of moment equations necessary to have a closed system depends on the kinetics of nucleation and crystal growth. For studying the properties and behaviour of the system a second order moment equation reduction is used in simulation. 15 Simulation results and discussion Applying the power law forms of kinetic expressions, we assume that the exponents are: g=1, b=2 and j=1. Then ( ) ( ) ⎟⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ −⎟ ⎠ ⎞ ⎜ ⎝ ⎛= − = = ∑ kks k k kb b ckt k tB 2,0,0 2 0 1 1 2)( πμ π 0 (56) i.e. nucleation describes also the possible dependence on the crystalline magma. As regards the crystal growth we have ( ) ( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ −⎟ ⎠ ⎞ ⎜ ⎝ ⎛= − = = ∑ jjs j j jg cj gktG 1,0,0 1 00 1 π π (57) while the expressions are give nin Eq. (55). Since the support of random parameters ωk is the compact interval [0,1] the micromixing parameters is characterised by the beta distribution ⎪ ⎩ ⎪ ⎨ ⎧ ≥≤ <<− = −− 1ifand0if,0 10if,)1( ),( 1 )( 11 kk k q k p k k qpf k ωω ωωω βωω k = a, b, c (58) including the limiting degenerate distributions given by Dirac-delta density functions as well. Note that, since we focus our attention on the effects of micromixing, in the model (47) crystals and agglomerates are treated as similar particles what is, in principle, a simplification. The mixing state of the feed consisting of species A and B can also be described by the beta distribution. We assume that the feed contains also species C having uniform concentration cc, in = cs where cs denotes the solubility of C. If species A and B in the feed are totally segregated with concentrations ca, in and cb, in and volumetric ratio φ, then the inlet population density function takes the form ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )inccinbbain inccbinaaincbain ccccc ccccctcccp ,,, ,,, 1 ,,, −−−+ +−−= δδδπφ δδδφπ 0 0 (59) so that the inlet joint moments are given as π1, 0, 0, in = π0, inφca, in, π0, 1, 0, in = (1 – φ)π0, incb, in π0, 0, 1, in = π0, incc, in, 2,,,0,0,2 inainin cφππ 0= ( ) 2,,,0,2,0 1 inbinin cφππ −= 0 , 2,,,2,0,0 incinin c0ππ = π1, 1, 0, in = 0, π1, 0, 1, in = π0, inφca, in cc, in π0, 1, 1, in = (1 – φ)π0, incb, in cc, in. The system of ODEs of the second order model equation reduction was solved in MATLAB environment using the process and kinetic parameters presented in Table 2. Simulation runs have been carried out to investigate the influence of micromixing and the interactions of micromixing and chemical reaction on the properties of crystal population produced. Table 2: Process and kinetic parameters used in simulation of reaction crystallizer 10=t s kb=1.0·10 8 Scol=10 s -1 kg=5.0·10 -11 cAin=2.2 mol dm -3 ρ=1500 g dm-3 cBin=2.2 mol dm -3 g=1.0 cCin=0.5 mol dm -3 b=2.0 cs=0.1mol dm -3 j=1.0 Sa=10 6 s-1 kr=4.0 10 1 mol dm-3 s-1 vf=1.0·10 -15 vn=1.0·10 -15 Effects of micromixing on the number of particles and the total mass produced in the crystallizer for three different reaction rates are presented in Figs 7 and 8. Here, the mixing state on microlevel is characterised by the expected values of parameters ωa, ωb and ωc defined as ∫ == 1 0 ,,,)( cbakdfm kkk kk ωωω ωω (60) Increasing rate of mixing on microlevel induces increased number of particles, crystals and aggregates, but at the same time the amount of solid mass is also increased. Under such circumstances the mean size of particles becomes 10-6-10-5 m of orders of magnitude. The diagrams in Figs 7 and 8 indicate also that the sensitivity of variation is increased with increasing rate of the chemical reaction. Indeed, these trends are justified by Figs 9 and 10 showing variations of the number of particles, as well as of the total volume of solid phase as a function of the rate of chemical reaction at different levels of micromixing. In these simulation runs the expected values of the micromixing parameters were constant mωa = mωb = mωc = 0.5. Figure 7: Influence of the expected value of micromixing parameters on the number of particles produced kr=1.0 kr=40 kr=0.4 μ 0 ( ) cba mmm ωωω == 16 Figure 8: Influence of the expected value of micromixing parameters on the mass of particles produced Figure 9: Influence of the rate of chemical reaction on the number of particles produced Figs 9 and 10 show that interactions between the chemical reaction and micromixing appear to be significant in the interval kr ∈ [10 -2, 102]. When kr > 10 2 then micromixing is the rate limiting process while in the case of kr < 10 2 the chemical reaction turns to be the rate limiting process. Conclusions The population balance approach provides a good framework for modelling crystallization systems. In extended form [12, 15, 20], it takes into consideration not only changes of the particles themselves, i.e. nucleation, growth, breakage and aggregation, but also changes of different entities carried by the particles, i.e. mass of chemical species and heat, induced by direct interparticle interactions. This extension often requires formulating multi-dimensional population balance models, and allows to model also interactive multi- populations involving different kinds of disperse elements. Figure 10: Influence of the rate of chemical reaction on the mass of particles produced The method of moments that seems to be an adventageous tool of handling one-dimensional population balance models can be applied also for multi-dimensional models defining the joint moments of variables describing simultaneously more internal properties of the dispersed elements. For the sake of illustration, a two-dimensional model was presented for describing a cooling crystallizer producing needle-shape crystals, and a two-population model was developed for describing micromixing of solution in reaction crystallisation. In the latter case, the population balance equation of fluid elements turned to be three-dimensional by itself. The second order moment equation reductions for joint moments of internal variables, which were closed by means of cumulant neglect closures, proved to be very useful in handling the rather complex models. NOTATION 〈c〉 mean (macrolevel) concentration, kg/m3 〈G〉 mean crystal growth rate, m/s 〈L〉 mean crystal size, m 〈B〉 mean nucleation rate, no/m3s t mean residence time, s a heat transfer surface in a unit volume of suspension ak coefficients of solubility (k = 0,1,2) b exponent of nucleation rate B nucleation rate, no/m3h c concentration, kg/m3, kmol/m3 Dab dimensionless parameter for secondary Dap dimensionless parameter for primary nucleation G crystal growth rate, m/s g exponent of crystal growth rate j exponent of nucleation rate kb coefficient of nucleation rate, (no/m 3s)⋅(m3/kg)-b kg coefficient of crystal growth rate, (m/s)⋅(m 3/kg)-g kr reaction rate coefficient, l/mol s kV volumetric shape factor L crystal size, m n population density function of crystals nucleation kr=40 kr=0.4 kr=1.0 μ 0 μ 1 μ 1 ( ) cba mmm ωωω == 5.0=== cba mmm ωωω 5.0=== cba mmm ωωω kr kr 17 p population density function of fluid elements R rate of chemical reaction, kmol/m3s s scale factor U heat transfer coefficient, W/m2K T temperature, C, K t time, s V volume, m3 vf volume of fluid elements, m 3 x dimensionless moments y dimensionless concentration GREEK LETTERS α dimensionless parameter ρ density of crystals [kgm-3] θ dimensionless nucleation rate ξ dimensionless time ε void fraction μm,k (m,k) th joint moment of crystal sizes πm,k m,k) th joint moment of concntrations INDICES a,b,c chemical species A,B,C b secondary nucleation in input out output p primary nucleation s equilibrium saturation concentration REFERENCES 1. HULBURT H. M., KATZ S.: Chemical Engineering Science, 1964, 19, 555-574 2. RANDOLPH A. D., LARSON M. A.: Theory of Particulate Processes (Second Ed.), Academic Press, New York, 1988 3. KAFAROV V. V., DOROKHOV I. N, KOLCOVA E. M.: Sistemnyj analiz processov khimicheskoj technologii. Processy massovoj kristallizacji iz rastvorov i gazovoj fazy. Nauka, Moskva, 1983 4. TAVARE N. S.: Industrial Crystallization. Process Simulation, Analysis and Design. Plenum Press, New York, 1995 5. BLICKLE T., LAKATOS B. G, MIHÁLYKÓ CS.: Szemcsés rendszerek matematikai modelljei. A kémia újabb eredményei 89. Akadémiai Kiadó, Budapest, 2001 6. HIDY G. M., BROCK J. R.: The Dynamics of Aero- colloidal Systems. Pergamon Press, New York, 1970 7. PANDIS S. N., SEINFELD J. H.: Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, Wiley, New York, 1998 8. FRIEDLANDER S. K.: Smoke, Dust, and Haze: Fundamentals of Aerosol Dynamics (2nd Ed.), Oxford University Press, Oxford, 2000 9. GIDASPOW D.: Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions. Academic Press, Boston, 1994 10. RAMKRISHNA D.: Population Balances. Theory and Applications to Particulate Systems in Engineering. Academic Press, San Diego, 2000 11. LAKATOS B. G., MIHÁLYKÓ CS., BLICKLE T.: Modelling of interactive populations of disperse systems. Proc. 2nd International Confrerence on Population Balance Modelling, Valencia, 2004, pp.72-76 12. LAKATOS B. G., MIHÁLYKÓ CS., BLICKLE T.: Chemical Engineering Science, 2006, 61, 54-62 13. MIHÁLYKÓ CS., LAKATOS B. G., MATEJDESZ A., BLICKLE T.: International Journal of Heat and Mass Transfer, 2004, 47, 1325-1334. 14. SÜLE Z., MIHÁLYKÓ CS., LAKATOS B. G.: Modelling of heat transfer processes in particulate systems. Computer-Aided Chemical Engineering 21A, Elsevier, Amsterdam, 2006, pp. 589-594 15. LAKATOS B. G., SÜLE Z., MIHÁLYKÓ CS.: International Journal of Heat and Mass Transfer 2008, 51, 1633-1645 16. SÜLE Z., MIHÁLYKÓ CS., LAKATOS B. G.: Chemical and Process Engineering, 2008, 29, 201-213 17. LAKATOS B. G., ULBERT ZS., BLICKLE T.: Population balance model for micromixing in continuous crystallizers. Proc. 14th International Symposium of Industrial Crystallization. Cambridge, UK, 1999 18. LAKATOS B. G., ULBERT ZS.: Population balance model for micromixing in MSMPR reaction crystallizers. Proceedings 3rd ECCE, Dechema, Nuremberg, Germany, Paper 374, 2001 19. ULBERT ZS., LAKATOS B. G.: Chemical Engineering Science, 2005, 60, 3525-3538 20. LAKATOS B. G.: Chemical Engineering Science, 2008, 63, 404-423 21. POHORECKI R., BALDYGA J.: Chemical Engineering Science, 1983, 38, 79-83 22. TAVARE N. S.: Chemical Engineering Science, 1994, 49, 5193-5201 23. TAVARE N. S.: AIChE Journal, 1995, 41, 2537-2548 24. FALK L., SCHAER E.: Chemical Engineering Science, 2001, 56, 2445-2457 25. LAKATOS B. G., NAGY G.: Probability function based micromixing model of crystallizers. Industrial Crystallization’87. (Ed. by Nyvlt J., Zacek S.), Praha: Academia, 325-328, 1987 26. ULBERT ZS., PhD thesis. University of Veszprém, Veszprém, Hungary, 2003 27. MA H., TAFTI D., BRAATZ R.: Industrial and Engineering Chemistry Research, 2002, 41, 6217- 6223 28. BRIESEN H., Chemical Engineering Science, 2006, 61, 104-112