CRYSTAL GROWTH & DESIGN
Semiempirical Model of a Batch Crystallizer
2003 VOL. 3, NO. 5 733-739
M. Parisi* and A. Chianese Department of Chemical Engineering, University of Rome ,La Sapienza., via Eudossiana 18, I-00184 Rome, Italy Received August 2, 2002
ABSTRACT: An approach to simulate the operation of batch crystallization of brittle materials by means of a semiempirical model is proposed. The model is derived from an experimental plan specifically designed to reduce the correlation among the values of the unknown parameters to be determined. The experimentation was carried out by using a pilot-scale DT crystallizer, 13.5 L in capacity, with and without seeding. The examined system was the crystallization of potassium sulfate from aqueous solutions. The developed model, accounting for the first nucleation point, crystal attrition, and growth rate dispersion, gave satisfactory results. The work, thus, indicates an overall procedure to provide a prediction model for an existing batch crystallizer. Introduction Batch cooling crystallizers are widely used at the industrial scale. The crystal size distribution (CSD) of the solid product is determined by several factors: initial nucleation point, hydrodynamics of the suspension, run-time, nucleation and crystal growth rates. When seeding is adopted, the mass and size of the seeds also affect the final CSD. As a consequence, modeling a batch crystallization operation is a very hard task. An industrial crystallizer is usually designed by scaling up or down other existing equipment operating the same crystallization process. In other words, a simulation model is very seldom adopted by the crystallization equipment suppliers for the design of a new crystallizer because of the great complexity of a physical prediction model. The weakest point of such a prediction model is the estimation of the nucleation point and of the secondary nucleation rate. All the nucleation phenomena are in fact very dependent on the hydrodynamics and on the geometry of the specific crystallizer. Gahn and Mersmann1 have recently proposed quite a sophisticated model to predict the generation rate of fragments on the basis of the geometry of a DT crystallizer and the impeller rotation speed. As a matter of fact, Bermingham et al.2 have recently noted the difficulty of simulating an ammonium sulfate pilot-scale crystallizer by means of this approach. For this reason, it seems, at the moment, more realistic to develop crystallizers simulation models that are derived and validated from experimental data more than predicting the crystallizer performances on the basis of general knowledge of the single crystallization phenomena. Along these lines are the recent approaches of adopting empirical models, which apply neural network techniques to evaluate the performances of an existing crystallizer. Giulietti et al.3 adopted a batch propagation neural network to simulate a batch crystallization for the production of copper and zinc sulfate from aqueous solution. However, a fully empirical model, such as the one mentioned above, is not robust enough when applied outside the parametric space adopted for its training. In this work, a semiempirical model of a batch crystallizer based on a physical description of single crystallization phenomena is presented. To determine
Figure 1. The experimental apparatus.
the kinetic coefficients, the design of a suitable experimental plan is proposed: the adopted approach might be applied during the test run of an industrial batch crystallizer, since it leads to a satisfactory description of the crystallization operation on the basis of a low number of runs to be performed. The model accounts for all the main phenomena occurring in a mixed suspension crystallizer. In particular, the generation of fragments due to the collisions of brittle crystals and the growth rate dispersion is considered. The batch crystallization of potassium sulfate from aqueous solution in a pilot-scale equipment is examined. Experimental Section Apparatus and Procedures. The adopted experimental apparatus is shown in Figure 1. It consists of a 13.5 L capacity
10.1021/cg0255669 CCC: $25.00 © 2003 American Chemical Society Published on Web 07/19/2003
734 Crystal Growth & Design, Vol. 3, No. 5, 2003
Parisi and Chianese
Table 1. Operating Condition of the Experimental Runs
unseeded runs standard seeded runs overseeded runs
ID
T0, °C
dT/dt, °C/h
w0, kg of salt/ kg of water
A B C D E F
42 56 42 56 42 42
7 14 7 14 7 7
0.155 0.177 0.153 0.175 0.153 0.153
seeding kg of salt/kg of water
9.93 × 10-4 3.70 × 10-3 1.11 × 10-2 3.33 × 10-2
size, µm
NO NO 300-355 300-355 850-1000 850-1000
Table 2. Main Experimental Results ID
final mass of crystals, kg of salt/kg of water
wf, kg of salt/ kg of water
Lwm, µm
fines % (L < 255 µm)
solute unbalance %
A B C D E F
2.10 × 10-2 4.81 × 10-2 1.96 × 10-2 5.03 × 10-2 2.23 × 10-2 5.75 × 10-2
0.135 0.128 0.136 0.131 0.133 0.131
454 536 756 862 1240 1140
15.4 12.8 5.4 3.3 3.6 1.1
-0.7 -0.1 -1.0 -1.2 -1.0 -1.4
batch crystallizer with a draft tube (D) and four vertical baffles (B); the bottom of the crystallizer was roundly shaped to avoid the segregation of crystals. Agitation was provided by means of a tree blade marine propeller (P) rotating at a constant speed of 940 rpm; the calculated dispersed power was around 0.4 W/kg. The temperature inside the crystallizer was continuously detected and manually controlled by means of a cooling coil (C) placed inside the draft tube and connected to a thermostatic bath. The adopted experimental procedure was as follows. A saturated solution of potassium sulfate was prepared at the initial temperature of the run, 42 or 56 °C, and then the solution was cooled at a constant cooling rate of 7 or 14 °C/h. All the runs had a duration time of 2 h and a final temperature of 28 °C. Seeding was adopted in some of the runs. The seeds, produced in previous batch runs, were added at a temperature just lower than the saturation one. Throughout each run, every 20 min a sample of the mother liquor was withdrawn and the salt concentration was measured by dryness. The sampling was carried out by means of a syringe fitted with a filter at its entrance. At the end of each run, the crystals were separated from the mother liquor by filtration and then air-dried. The crystal size distribution was determined up to 106 µm by means of a laser light scattering based instrument and by sieving for larger dimensions. Three different kinds of runs were carried out (see Table 1): (i) Runs without seeding (runs A and B). (ii) Runs with a standard amount of seeds in the size range 300-355 µm (runs C and D). (iii) Runs with a greater amount of seeds in the size range 850-1000 µm (runs E and F). Seeds of a larger size were used in the fragmentation runs to give evidence of the secondary nucleation by collision mechanisms. Potassium sulfate at a purity grade higher than 99.5% b.w. was used. Experimental Results. The main performances of the batch crystallization runs are reported in Table 2, whereas the supersaturation profiles are plotted in Figure 2a,b. The effects of the operating conditions were as expected. In particular: (i) The runs carried out without seeding (A and B) exhibited the lowest weight mean crystal sizes and the highest percentage of fines. (ii) The addition of seeds drastically decreased the amount of fines from more than 12% (runs A and B) down to less than 6% (runs C, D, E, and F). (iii) The runs operating with a wider temperature cooling range (runs B and D) produced a higher mass of crystals.
Figure 2. Supersaturation vs time of the experimental runs: (a) runs A, B, C, and D; (b) runs E and F. (iv) The run F, performed with the highest mass of seeds, gave rise to the lowest percentage of fines, due to the small supersaturation attained throughout the run, less than or equal to 0.06 (see Figure 2b). The supersaturation profiles are shown in Figure 2. They were calculated with respect to the equilibrium data of the potassium sulfate in water solution, that are known from the literature4 and can be expressed by the following relations:
weq ) 7.4 × 10-2 + 1.88 × 10-3T
for T < 43.8 °C
weq ) 8.8 × 10-2 + 1.56 × 10-3T
for T g 43.8 °C
(1)
The supersaturation profiles are very different for the various runs. Obviously, the maximum supersaturation is achieved by the runs carried out without seeding. For these runs, the abscissa of the peak in the supersaturation profile approximately corresponds to the nucleation point, which is
Semiempirical Model of a Batch Crystallizer
Crystal Growth & Design, Vol. 3, No. 5, 2003 735
approached earlier when the cooling rate is lower. For the seeding runs, the supersaturation values decreases on increasing the seed mass as expected. When the seed mass is equal to or less than 50 g the presence of a peak in the supersaturation profile is evident, whereas at higher values a peak cannot be clearly identified in the supersaturation profile, which is quite flat in the second half of each run.
Simulation Model The modeling of a batch crystallizer is based on three equations representing the unsteady state heat, mass, and population balances.5 When the crystallizer’s temperature is strictly controlled, the heat balance can be substituted by the temperature profile, which in case of constant cooling rate is as follows: T(t) ) T0 - Rct
(2)
where T is the temperature in the crystallizer, t is the time, Rc is the adopted constant cooling rate, -dT/dt, and the subscript 0 refers to the initial conditions. The solute mass balance can be expressed by the following equation: MT0 + w0 ) MT(t) + w(t)
(3)
where MT(t) and w(t) are the magma density and the solute concentration, respectively, both expressed as kg of solute/kg of solvent. The population mass balance, according to Randolph and Larson,6 can be expressed by the equation: δ(Gn) δn )+ b(L) - da(L) δt δL
(4)
where L is a characteristic dimension of the crystals, G is the linear growth rate equal to dL/dt, n(L) is the population density function, and b(L) and da(L) are the birth and death density functions, respectively. Equation 4 can be discretized and numerically integrated with the boundary condition at t ) 0: n(L) ) ns(L) n(L) ) 0
for seeded runs
for unseeded runs
(5)
where ns(L) is the population density function of the seeds. The birth rate of nuclei is determined by the primary heterogeneous nucleation or secondary nucleation by catalytic mechanism, respectively, in the absence and presence of crystals in suspension. Primary nucleation occurs when the undercooling overcomes a threshold value that depends on the system and on the cooling rate. When crystals are present in suspension, the value of the threshold undercooling that allows the formation of new nuclei is reduced due to the catalytic effect of the crystals surface.4 The nucleation rate was considered equal to zero before the first nucleation occurred. The following empirical expression was adopted to determine both primary heterogeneous and secondary catalytic nucleation:7 B(L f 0) ) kn(1 + RMT)σ i for t g tn
(6)
where kn, R, and i are empirical parameters and tn is the nucleation point corresponding to the time when first nucleation occurs. In absence of crystals MT is equal to zero and primary heterogeneous nucleation occurs. Furthermore, fragments are generated due to secondary nucleation by crystal collision mechanisms. According to Chianese et al.,8 the attrition and abrasion of brittle crystals can produce fragments of up to about 150 µm. The size distribution of the generated fragments can be expressed by means of the number density function, defined as the ratio between the fraction of fragments in the size range ∆L and the range width, for ∆L f 0: fL(L) )
1 dNFT(L) NFT dL
(7)
where NFT is the total number of fragments generated per unit volume and per unit time in the suspension. The fragment density function, fL(L), is well-represented by the following exponential equation: fL(L) ) ALb
(8)
where, for brittle material, the exponent b is equal to -2.96.8 The constant A can be determined by applying the normalization property of the density function between the minimum and maximum dimension of the fragments: A)
b+1 b+1 Lb+1 F,min - LF,max
(9)
Thus, the birth density function for the fragments of size L is b(L) ) NFT fL(L)
(10)
The total number of fragments generated per unit volume and per unit time, NFT, can be expressed by the equation NFT )
∫
∞
0
fn(L)n(L) dL
(11)
where fn(L) is the average number of fragments generated by collision from a crystal of size L per unit time. The function fn(L) depends on several factors, as crystal size and concentration of crystals in suspension, supersaturation, dispersed energy per unit volume, and per unit time, etc. In this work, a constant impeller speed was adopted, thus the dependence on the power dispersed by the impeller was neglected, and the following empirical expression was assumed: fn(L) ) kFkVL3(1 + Cσ 2)MT
(12)
where σ is the relative supersaturation, MT is the magma density, kF and C are empirical parameters and kV is the volume shape factor, defined by the equation: Vc ) kVL3
(13)
The volume shape factor of the potassium sulfate crystals was estimated as a function of the crystal size by the relationships proposed by Budz et al.9
736
Crystal Growth & Design, Vol. 3, No. 5, 2003
Parisi and Chianese
During the preliminary simulation work, the term Cσ 2 resulted much bigger than 1; thus, eq 12 was simplified as: fn(L) ) k′FkVL3σ 2MT
(14)
where k′F ) kFC. The death density function da(L) can be represented by means of a linear attrition rate, Ga ) -dL/dt, acting in the opposite direction with respect to that of the linear crystal growth rate G, as proposed by Polish et al.:10 -da(L) )
δ(nGa) δL
(15)
The linear attrition rate Ga(L) can be related to the function fn(L) and to the mean volume of the fragments VF,mean, defined as: VF,mean )
∫
LF,max
LF,min
kVL3fL(L) dL
(16)
dVc ) VF,mean fn(L) dt
dL VF,mean fn(L) ) dt 3k L2
6 6 6 13 20 20 21 22
1.4 × 0.4 6 1.0 × 108 5 50 110 0.2
1010
unit s-1
no. m-1 dimensionless dimensionless no. s-1 m-3 dimensionless dimensionless m s-1 dimensionless
having a growth rate between G - ∆G/2 and G + ∆G/2 is f(G/Gmean) )
∆G/Gmean(k - 2)k-1e-(k-2)/(G/Gmean) (G/Gmean)kΓ(k - 1)
(20)
where k is the parameter of the inverse gamma function and Gmean is the mean growth rate expressed by a second-order growth rate equation: Gmean ) k′Gσ 2e-Ea/RT
(21)
(18)
GF,mean ) rgGmean
the linear attrition rate can be determined by the equation: Ga(L) ) -
value
(17)
and expressing the volumetric attrition rate as: -
kn a i k′F k (new crystals) k (seeds) k′G rg
equation
where k′G is the growth rate coefficient and Ea is the activation energy, equal to 40 400 J/mol for the potassium sulfate.14 The growth rate order equal to 2 is the usual value adopted for potassium sulfate crystallization.7,14 Many authors reported that the fragment mean growth rate is considerably lower than that of the nuclei.15 Accordingly, the following expression for the mean growth of fragment was assumed:
In fact, differentiating both the members of eq 13 with respect to time, we have dL 1 dVc ) dt k 3L2 dt V
Table 3. Parameters of the Models parameter
(19)
V
One of the main physical phenomenon to be taken into account in the simulation of a crystallization operation is the growth rate dispersion:11 crystals of the same material grow at different rates regardless of the same operative supersaturation and temperature. This phenomenon has been ascertained by the experimental work of Sherwood et al.12 Two different models are proposed in the literature to interpret the growth rate dispersion: the random fluctuation model,11 based on the hypothesis that the growth rate of all the crystals randomly fluctuates around a common mean value, and the constant crystal growth model, originally proposed by Larson et al.13 According to this latter model, the growth rate of a single crystal is univocally determined by the applied operating condition, but it may differ from crystal to crystal. This model is based on the observation that the growth rate of each crystal is related to the density of the dislocations over its surfaces and is constant throughout its life. In this work, the constant crystal growth mechanism was assumed for the growth rate dispersion. The function adopted to express the crystal growth distribution was the inverse gamma function.13 Assuming a discretized form of this distribution, the fraction of crystals
(22)
where the parameter rg, less than 1, was evaluated from the best fitting of the experimental runs. Numerical Solution of the Simulation Model. The simulation model was numerically solved by writing a computer code in Visual Basic. The crystal growth dispersion was taken into account by considering different classes of crystals, each having a constant crystal growth rate.17 Growth rate classes were considered and for each class a specific population balance was written. Different population balance equations were considered for the fragments and seeds. In fact, for the fragments a different growth rate was adopted as mentioned above, whereas for the seeds a different lower inverse gamma distribution was adopted for the growth rate dispersion. The reason is that the seeds are crystals produced in a narrow size range, which thus belonged to a particular fraction of the crystal population with a narrow growth rate dispersion. Each partial differential population equation was solved by a finite difference procedure, by adopting a crystal size range of 2 microns and an increasing time integration interval. For the integral-differential set of equations, a stable solution was obtained. The computation time by using a Pentium III GHz computer was between 5 and 20 s. Results and Discussion The adjusting parameters of the model, reported in Table 3, were obtained by the best fitting of the CSDs curves.
Semiempirical Model of a Batch Crystallizer
Figure 3. Mass fraction distribution for run A.
The following procedure was applied to estimate a robust and quite approximate initial set of parameters: (i) The nucleation rate and the growth rate coefficients were evaluated on the basis of the measured supersaturation profiles. First of all, from the maximum of the supersaturation profiles of the two unseeded runs the values of the metastable zone width ∆TMAX were derived and attributed to the two adopted cooling rates, dT/dt. By applying the Ny´vlt method8 on the couples of values (∆TMAX, dT/dt) the nucleation parameters kN and I were estimated. Then, the third nucleation parameter R, accounting for the increase of the nucleation rate due to the magma density of crystals, was determined by merging the percentage of fines obtained from the seeded runs C and D. Finally, the growth rate coefficient k′G was calculated by the displacement of the median size of the seeded crystals in the runs C and D. (ii) The two growth rate dispersion coefficients referring the new generated crystals and seeds, respectively, were determined by comparing the experimental and calculated size distributions of the two classes of crystals in the runs A, B, C, and D. (iii) The generation rate and the growth rate reduction factor of the fragments, i.e., kF and rg, respectively, were estimated by fitting the CSDs for the runs E and F, operated at a high seed mass. The first approximation values of the parameters were then tuned by the best fitting of CSDs, in terms of cumulative mass density, for all the runs. It is important to notice that the parameter values were not strongly modified by the tuning procedure: in fact, the parameters concerning the growth rate determined by fitting the runs A, B, C, and D, were able to describe also the runs E and F. In Figures 3-5 the experimental values of CSDs are compared with the calculated profiles. A satisfactory agreement between the experimental values and those calculated by the model proposed in this work was obtained for all the runs. Figures 6-8 show in a semilogarithmic plot the comparison between experimental and predicted profiles of crystal population density distributions. The curves calculated by the model well fit the experimental data for all the runs. The only significant deviation between the predicted and the experimental values of the size population density distribution concerns the fine crystals for runs A and B. This happens
Crystal Growth & Design, Vol. 3, No. 5, 2003 737
Figure 4. Mass fraction distribution for run C.
Figure 5. Mass fraction distribution for run E.
Figure 6. Population density distribution for run A.
Figure 7. Population density distribution for run C.
also because the fitting of the experimental results did not refer to the population density, but to the mass fraction distribution. In this work, the mass distributions instead of the logarithm of the population density distribution was used to determine the adjusting parameters of the model because this choice enables the achievement of more
738
Crystal Growth & Design, Vol. 3, No. 5, 2003
Parisi and Chianese
a rather sophisticated model has to be adopted. The dispersion growth rate assumption is required to account for the size distribution of coarse crystals, whereas the distribution of fines are well predicted only if breakage phenomena and the distribution of the generated fragments are properly considered. A model, which takes into account all these phenomena as that developed in this work, is flexible, thus allowing the interpolation of monomodal and bimodal crystal size distributions and of supersaturation profiles exhibiting a sharp peak (unseeded runs) and a quite smooth trend (runs with high seed mass). Figure 8. Population density distribution for run E.
List of Symbols information from the experimental data. First of all, the use of the logarithmic scale instead of the linear scale reduces the deviations among the experimental data and get less sensitive the model response with respect to the parameters values. Moreover, the use of the population density distribution does not allow to clearly evaluate the wideness and the displacement of the peaks of the size distributions of the fines and coarse crystals. As a matter of fact, this latter information was very useful to evaluate the initial parameter set in this work. To make evidence of the importance of the growth rate dispersion a simulation of run C was made by adopting a constant value for the crystal growth equal to the mean growth rate. The simulation results reported in Figure 4 show a very narrow and high peak for the coarse crystals, strongly in contrast with the experimental size distribution. It has to be noted that the adopted model is capable of simulating quite well batch crystallization runs performed at very different operating conditions; thus, it is flexible enough to account for a number of phenomena, i.e., primary nucleation, attrition, growth rate dispersion, which occur to a different extent in the different runs. The adjusting parameters, reported in Table 3, are not physically meaningful, and for this reason the described model has to be considered a semiempirical one. When a wider experimental data set is available the parameters may be better tuned and the model would achieve more robustness. However, some parameters are in good agreement with those derived from a experimental investigations on single phenomena. This is the case for the growth rate function. Mydlarz and Jones16 report, for T ) 28 °C and σ ) 0.097, values of the growth rate of potassium sulfate, depending on the crystal size, in the range of 0.6-1.5 × 10-7 m/s, for crystals size between 100 and 1000 µm. For the same operating conditions, the calculation of the mean growth rate adopted in this work gives a value of 1.2 × 10-7 m/s, which is in very good agreement with the value reported by Mydlarz and Jones. Conclusion This work shows how it is possible to develop a simulation model of a batch crystallization process when an appropriate experimental plan is made. To evaluate with satisfactory accuracy the size distribution of the produced crystals from a seeded or unseeded batch run
b(L) da(L) fL(L) fn(L) G Ga i k′F kG kn kv L MT n NFT Rc rg t T V w
birth density function, no. s-1 m-1 death density function, no. s-1 m-1 number density function of the fragment, no. m-1 number of fragments generated by collision per unit volume and per unit time by a parent crystal of size L, no. m-3 s-1 growh rate, m s-1 attrition rate, m s-1 parameter of the birth density function, dimensionless parameter of the fragment generation function, no. s-1 m-3 parameter of the growth rate function, m s-1 parameter of the birth density function, no. s-1 m-1 volume shape factor, dimensionless linear dimension of the crystals, m magma density, kg of crystalline substance in suspension/kg of solute population density function, no. m-1 number of fragments generated per unit volume and per unit time, no. m-3 s-1 cooling rate, K s-1 fragment growth rate reduction factor, dimensionless time, s temperature, K volume of the particles, m3 solute concentration, kg of solute/kg of solvent
Greek Letters R parameter of the birth density function, dimensionless σ relative supersaturation, dimensionless Subscripts 0 initial c crystal eq equilibrium F fragment max maximum value mean mean value min minimum value wm weight mean size
References (1) Gahn, C.; Mersmann, A. Chem. Eng. Sci. 1999, 54, 12731282. (2) Bermingham, S. K.; Neumann, A. M.; Verheijen, P. J. T.; Kramer, H. J. M. BIWIC 2001, Delft, September 2001, pp 118-125. (3) Giulietti, M.; Terra, L. R.; Guardani, R. Proceedings of 14th Int. Symp. on Industrial Crystallization; Cambridge, 1999. (4) Mullin, J. W. Crystallization; Butterworth Heinemann: Oxford, 1993. (5) Parisi, M.; Chianese, A. BIWIC 2001, Delft, September 2001, pp 126-133.
Semiempirical Model of a Batch Crystallizer (6) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes; Academic Press: New York, 1988. (7) Ny´vlt, J.; So¨hnel, O.; Matuchova`, M.; Broul, M. The Kinetics of Industrial Crystallization; Elsevier: Oxford, 1985. (8) Chianese, A.; Sangl, R. G.; Mersmann, A. B. Chem. Eng. Commun. 1996, 146, 11-12. (9) Budz, J.; Jones, A. G.; Mullin, J. W. Ind. Eng. Chem. Res. 1987, 26, 820-824. (10) Polish, J.; Mersmann, A. Annual meeting of German Chemical Engineering in Stasbourg, 1986. (11) Randolph, A. D.; White, E. T. Chem. Eng. Sci. 1977, 32, 1067-1076.
Crystal Growth & Design, Vol. 3, No. 5, 2003 739 (12) Sherwood, J. N.; Ristic, R. I. Chem. Eng. Sci. 2001, 56, 2267-2280. (13) Larson, M. A.; White, E. T.; Ramanarayan, K. A.; Berlunglund, K. AIChE J. 1885, 31. (14) Jones, A. G.; Budz, J.; Mullin, J. W. AIChE J. 1886, 32, 2002-2008. (15) Bravi, M.; Chianese, A. ISIC 98, Tianjin, 1998, 210-215. (16) Mydlarz, J., Jones, A. G. Chem. Eng. Comm. 1990, 90, 47-56.
CG0255669
