6228
Ind. Eng. Chem. Res. 2001, 40, 6228-6235
Adaptive Mesh Method for the Simulation of Crystallization Processes Including Agglomeration and Breakage: the Potassium Sulfate System Gibaek Lee* Department of Industrial Chemistry, Chung-Ju National University, Chung-Ju, Chungbuk 380-702, Korea
En Sup Yoon School of Chemical Engineering, Seoul National University, Seoul 151-742, Korea
Young-Il Lim, Jean Marc Le Lann, Xuaˆ n-Mi Meyer, and Xavier Joulia Laboratoire de Ge´ nie Chimique (LGC, UMR-CNRS 5503), INPT-ENSIACET, 18 Chemin de la loge, F-31078 Toulouse Cedex 4, France
This study focuses on the adaptive mesh method for simulation of crystallization processes. Because the models of the processes are very complex, giving partial and ordinary differential, algebraic, and integral equations, it is important to get a numerical solution. The simulation is based on the method of characteristics to avoid the difficulties of numerical diffusion and stability that occur in other methods (e.g., finite difference methods for first-order derivatives) used to solve these equations. The simulation includes the kinetics of breakage as well as agglomeration. This study combines an adaptive mesh method with the method of characteristics in order to improve the accuracy and efficiency of the simulation. The proposed method can handle stiff nucleation problems of the characteristic method and conserve the crystal number and mass during mesh adaptation. The proposed method is validated through the simulation of a potassium sulfate crystallization. Introduction Crystallization is one of the most important unit operations and is widely used in the chemical industry for purification and separation.1 An important property in crystallization is the crystal size distribution (CSD), which can be modeled in models by the population balance equation (PBE). However, because the equation is a hyperbolic partial differential equation coupled with other ordinary differential and algebraic equations, it is very difficult to get a solution. The kinetics of agglomeration and breakage in crystallization are expressed in integral equations. Often, the PBE with agglomeration and breakage cannot be solved analytically, and therefore numerical solutions are required. Furthermore, the recent model-based control of crystallization processes needs the numerical solutions to be accurate and available in real time. It should be noted that Tarvare et al.2 have obtained analytical solutions of the population balance for cooling or evaporation crystallization under the limiting conditions. A general numerical method to solve integro-partial differential algebraic equations (IPDAE) systems does not exist, though several numerical techniques, e.g., the
method of weighted residuals, method of moments, orthogonal collocation, and collocation on finite element, have been proposed.3 Because of their numerical simplicity, several investigations have been based on discretization methods that are similar to finite difference formulas.4,5 Among these methods, a combination method of the characteristics with the method of discretization has an excellent quality, although its implementation is difficult.6 In this study, we adopted the combination method and modified it to enhance the solution quality. The method will be introduced after the description of the crystallization model. Solving a partial differential equation (PDE) whose solutions display steep moving fronts has the need for numerical adaptation of the mesh.7 Because many meshes are required to retain accuracy at the front when the solution changes very rapidly, it needs an excessive amount of computation. To increase its efficiency and to guarantee accuracy of the result, we used an adaptive mesh method to divide or eliminate meshes adaptively. Finally, we present some numerical examples in order to illustrate the proposed method. The Model
* To whom correspondence should be addressed. Tel: +82-43-841-5230. Fax: +82-43-841-5220. E-mail: glee@ gukwon.chungju.ac.kr.
The mathematical equations describing the dynamic behavior for crystallization are given below. The model
10.1021/ie010443r CCC: $20.00 © 2001 American Chemical Society Published on Web 11/30/2001
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001 6229
includes the classical formulation expressing the mass and energy conservation as functions of time.
Partial mole balance equation: d(yiUg) d(xiUl) d(Us) + + ) dt dt dt
∑Fkxik
(1)
many researchers have assumed that the agglomeration kernel is size-independent. However, this study does not use that assumption. Among various agglomeration kernels suggested in the literature, Smoluchowski’s kernel is selected for the potassium sulfate system.10 It is proportional to the particle size, u, and v of the agglomerating particles. Therefore, a large particle can aggregate more frequently than a small particle.
Global mole balance equation: dUg dUl dUs + + ) dt dt dt
A(u,v) ) Kagg(u1/3 + v1/3)3
∑
Fkck
(4) Breakage: division of a big crystal into smaller ones. Like agglomeration, breakage involves two contributions to the population balance.
Energy balance equation: d(HgUg) d(HlUl) d(HsUs) + + ) dt dt dt
∑FkHk + Q - QC (3)
Randolph and Larson8 developed the following population balance concept in the crystallization process.
dn d log V d(Gn) +n + ) B0δ(v - v0) + Ba - Da + dt dt dv Fknk B b - Db + (4) V
∑
It is used to properly model a population of crystals suspended in a continuous phase. The four mechanisms involved in this equation are the following: (1) Nucleation: creation of new crystals from any one, or a combination, of primary and secondary nucleation or by attrition. The nucleation rate is often expressed as a power law model. The nucleation rate constant, Kn, depends on a number of variables.
B0 ) Kn∆C
b
(5)
(2) Growth: increase of the crystal size due to solute deposition. The overall growth rate is expressed by the empirical equation. The growth rate constant, Kg, depends on variables such as temperature, hydrodynamics, and impurities.
G ) Kg∆CgLlkve-Eg/RT
(6)
(3) Agglomeration: creation of a bigger particle from two smaller particles. The agglomeration of two particles of size s and v - s into a particle of size v has two effects on the population.9
Birth density function of particle of size v from particles of size s and v - s: Ba )
∫vvA2 n(v - s)n(s) ds 0
(7)
Death density function of particle of size v by agglomeration with particle of size s: Da )
∫v∞An(v) n(s) ds 0
(9)
(2)
(8)
The formulation of the agglomeration kernel, A(v-s,s), is chosen to correspond to the particular mechanism of agglomeration. Though many theoretical and empirical formulations of agglomeration kernels are available,
Birth of particles of size v by breakage of particles of bigger size: Bb )
∫v∞b(v,s) Bn(s) dx
(10)
Death of particles of size v by breakage into smaller particles: Db ) Bn(v)
(11)
To model breakage, this study uses the following form of kernel, where the breakage rate is proportional to the particle volume.
B(v) ) Kbrev
(12)
The equation for the solid holdup from the crystal size distribution is also an integral equation in addition to the above.
Us ) Fs
∫0∞kv(v) n(v) v dv
(13)
In addition, there are a number of algebraic equations expressing the physical phenomena such as supersaturation and vapor-liquid equilibrium as follows. The total model is a system of IPDAE.
Supersaturation equation: ∆C ) C - Ceq
(14)
Vapor-liquid equilibrium equation: yi ) kixi
(15)
Summation equation on fractions:
∑xi ) ∑yi
(16)
Suspension volume equation: U l Us + )V Fl Fs
(17)
Total volume equation: Ug Ul Us + + ) VT Fg Fl Fs
(18)
Method of Characteristics. The study uses the concept of the combined method of characteristics with
6230
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001
and mass, η (in eq 20) due to agglomeration is given as the following equation.
{
discretization techniques, which Kumar and Ramkrishna proposed.6 In their approach, the entire size range is partitioned into small size ranges as with the other discretization methods. The size interval between vi and vi+1 is called the ith mesh (Figure 1). The particle population in this size range is represented by a size si (called the representative size11), such that vi < si < vi+1. Therefore, the population of the ith mesh, Ni, is defined as follows.
Ni(t) )
∫v
i
n(v,t) dv
, si(t) e v e si+1(t)
s (t) - si(t) η ) i+1 v - si-1(t)
Figure 1. Mesh used in the method of characteristics.
vi+1
si+1(t) - v
(19)
Through the mathematical transformation procedure known as the “method of characteristics”, Kumar and Ramkrishna obtained the following equation from eq 4, using the Leibnitz rule.
(23) , si-1(t) e v e si(t)
si(t) - si-1(t)
v ) sj(t) + sk(t)
Similarly, ni,k, being the contribution to the population at the ith representative size due to the breakage of a particle of size sk, is given by the simple form
ni,k )
∫ss
i
v - si-1 b(v,sk) dv + si - si-1
i-1
∫ss
si+1 - v b(v,sk) dv si+1 - si (24)
i+1
i
This study assumes binary breakage. Because each broken particle forms two smaller particles, the breakage function is constrained by the following equation:
∫0s b(v,sk) dv ) 2 k
dNi(t) dt
M
)
∑
ni,kBkNk(t) k)1 igjgk
∑ j,k
- BiNi(t) +
(
Here, b(v,sk) is the discretized number fraction of particles broken from size sk into size interval v. Hill and Ng10,12 assumed that the size of the first mesh is zero and suggested the following breakage function.
)
1 1 - δj,k ηAj,kNj(t) Nk(t) 2
si-1e(sj+sk)esi+1
M
Ni(t)
∑ Ai,kNk(t) + ∫v k)1
vi+1 i
S(v) dv (20)
The total population of the ith size range is represented at a representative size si. The population density function can therefore be written as M
n(v,t) )
Ni(t) δ(v-si) ∑ i)1
(21)
b(v,sk) ) 2/sk However, the method of characteristics suggested by Kumar and Ramkrishna6 does not assume that the size of the first mesh is equal to zero. Therefore, the constraint from eq 25 cannot be satisfied and should be modified as follows:
∫ss b(v,sk) dv ) 2 k
1
The particle size (for the boundary size as well as the representative size) changes with time.
dsi ) G(si) dt
(25)
Here, s1 is the representative size of the smallest mesh. Therefore, the breakage function is as follows:
b(v,sk) ) (22)
Kumar and Ramkrishna proposed a strategy to conserve two moments of number and mass when a new particle is formed by breakage or agglomeration. If a newly formed particle does not have its size corresponding exactly to a representative size, moments cannot be conserved. Therefore, the particle is assigned to the adjoining representative sizes such that two moments of interest are exactly conserved. To preserve number
(26)
2 sk - 2s1
(27)
Also, the breakage function is constrained by the following symmetry requirement:
b(v,sk) ) b(sk-v,sk)
(28)
To satisfy eq 28, the particle of size sk should be broken into the interval of [s1, sk - s1] and ni,k should be calculated by the following equations.
{
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001 6231
s -v
s -s
∫ss s22 - s1b(v,sk) dv ) sk2- 2s11,
ni,k )
2
1
i)1
si+1 - si-1 v - si-1 si+1 si+1 - v b(v,sk) dv + s b(v,sk) dv ) , i ) 2 to k - 2 i si - si-1 si+1 - si sk - 2s1 sk - sk-2 s12 sk-1 v - sk-2 sk-s0 sk - v b(v,s ) dv + b(v,s ) dv ) , i)k-1 k k sk-2 s sk-1 s - s sk - 2s1 (sk - 2s1)(sk - sk-1) k-1 - sk-2 k k-1 sk - sk-1 s12 2s1 sk-s1 v - sk-1 b(v,s ) dv ) + , i)k k sk-1 s - s sk - 2s1 (sk - 2s1)(sk - sk-1) sk - 2s1 k k-1
∫ss
i
i-1
∫
∫
(29)
∫
∫
Though the approach overcomes the diffusion and stability problems of the earlier techniques, a new difficulty arises with the nucleation. Because the size moves with time, nucleation cannot be represented when nuclei are smaller than the smallest particle size of the whole range. To surmount this problem, Kumar and Ramkrishna used the strategy of adding a new mesh with zero particles whenever the lower boundary of the first mesh coincides with the nuclei size. With this method, nuclei number can be represented correctly, but mass may not be conserved. The method results in a number of very fine meshes due to newly added meshes, and the system size (or mesh numbers) increases with time. To resolve the problem, they used a technique to eliminate the fine grids. However, they did not propose any criteria, which will select the fine grids to be eliminated. To answer this question, we propose the strategy to combine this technique with the adaptive mesh method.
on the first and second derivatives of the solution can be used. They are similar to the arc length and curvature monitor functions for the moving mesh method. The arc length function makes the nodes concentrate on regions of steep slope and the curvature function on regions of the corner in the solution. We used the mesh function based on the first derivative of the solution as follows:
f(v,t) ) x(∆v)2 + (∆vNv)2 ) ∆vx1 + Nv2 ) ∆v
x1 + (∂N∂v )
2
(30)
The meshes will be adapted to equidistribute the mesh function. Therefore, the meshes having mesh function values larger than the following tolerance will be refined. M
Adaptive Mesh Method Traditionally, adaptive mesh methods have been used to improve the accuracy and efficiency of numerical approximations to PDE systems. Large gradients or discontinuities can occur, for example, in shock waves or moving fronts. To approximate the solution accurately in these ranges, it is often necessary to generate a fine mesh where the solution is changing steeply. Also, for reasons of efficiency, the mesh should be coarse where the solution is smooth. Adaptation can seem costly, but without any refinement, numerical calculations would be wasteful, or even worse, not resolve some important aspects of the solution satisfactorily.13 As the solution changes, the mesh must also change to adaptively refine ranges where the solution is developing sharp gradients and to remove unneeded points from ranges where the solution is becoming smoother. Thus, the mesh must have a dynamic behavior in much the same way as the solution. For the moving mesh method, the number of meshes is fixed. The computational coordinate transforms the nonuniform physical mesh into a uniform computational mesh so that the meshes equidistribute or approximately equidistribute using given monitor functions such as the arc length or curvature of the solution, and the meshes concentrate in regions of rapid variation of the solution.14 To answer the nucleation problem, we adopted an adaptive mesh method that is combined with the strategy of the moving mesh method. For the decision tolerance to adapt meshes, the mesh functions based
TOL )
∆vix1 + Nv,i2 ∑ i)1
(31)
M
Also, we can choose the relative mesh function in order to consider the size range where density functions are very small. M
TOLrel )
∆vi ∑ i)1
x( ) ( ) 1
vref
M
2
+
Nv,i
Nref
2
(32)
When the mesh function of the solution is much smaller than the tolerance, the meshes will be eliminated to reduce the calculation time. If the tolerance for mesh elimination is the same as the mesh division, mesh division and elimination will be performed so frequently that the integration may be unstable.7 Practically, the tolerance of mesh elimination should be larger than a tenth of the tolerance for mesh division in order to prevent this difficulty. To solve the problem of the nucleation mesh, Kumar and Ramkrishna6 introduced a new mesh with zero population before the first mesh at every time step. The representative size of the mesh has the nuclei size, and the smallest and largest sizes of the mesh are the same as the representative size (v0 ) v1 ) s1). The largest size and representative size move with time, but the smallest size does not move. Because the representative size of
6232
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001
Figure 3. Temperature profiles for natural and constant cooling scenarios.
Figure 4. Comparison of the experimental16 and simulated CSD for the natural cooling case without agglomeration and breakage. Figure 2. Graphical representation for (a) mesh elimination and (b) mesh addition.
the nuclei mesh differs from the nuclei size, the mass in the mesh may not be correctly represented. Therefore, the adaptation time step should be very short where the mass in the nuclei mesh is large. For this reason, this study adopted the strategy that keeps the mass ratio in the nuclei mesh to all meshes under a prespecified constant. That is, a new mesh is added when the mass in the nuclei mesh is no longer negligible compared to the total mass. The mesh adaptation is performed when a new mesh is added at nuclei size or when 10% of the total meshes require the adaptation. The resultant differential algebraic equation system method was compiled in Visual Fortran 6.5 and solved in the DAE solver DISCo.15 In the solver, the routine to determine the mesh adaptation is performed as a state event. One mesh has seven variables in the system (population number, size, volume shape factor, and birth/death density functions by agglomeration/breakage), so seven variables per mesh are added to or removed from the integration system when the meshes are adapted. Among these variables, the volume shape factor can be directly calculated from other state variables, and the population number is calculated to conserve the two moments corresponding to number and mass. Consider the division of the ith mesh as shown in Figure 2a. Before the ith mesh, one mesh is added to become the new ith mesh. As meshes are renumbered,
the old ith mesh becomes the new (i + 1)th mesh, and so on. (Figure 2a). For the conservation of number and mass, the populations of new ith and (i + 1)th meshes should satisfy the following equations to calculate number and mass (see eq 13).
Niold ) Ninew + Ni+1new
(33)
kv(siold)sioldNiold ) kv(sinew)sinewNinew + kv(si+1new)si+1newNi+1new (34) From these equations, the following equations are obtained.
Ninew )
new
Ni+1
kv(si+1new)si+1new - kv(siold)siold kv(si+1new)si+1new - kv(sinew)sinew
)
Niold
kv(siold)siold - kv(sinew)sinew kv(si+1new)si+1new - kv(sinew)sinew
(35) Niold (36)
Consider the elimination of the ith mesh as shown in Figure 2b. By mesh renumbering, the old (i + 1)th mesh becomes the new ith mesh, and the process continues (Figure 2b). Similar to the above case, the new populations of (i - 1)th mesh and the population of the new
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001 6233
Figure 5. Comparison of the experimental16 and optimized CSD considering agglomeration and breakage.
ith mesh satisfy the following equations to calculate number and mass.
Ni-1old + Niold + Ni+1old ) Ni-1new + Ninew (37) kv(si-1old)si-1oldNi-1old + kv(siold)sioldNiold + kv(si+1old)si+1oldNi+1old ) kv(si-1new)si-1newNi-1new + new
kv(si
)si
new
new
Ni
(38)
From these equations, we can obtain the following equations to calculate the new population numbers. i+1
Ni-1new )
Nkold[kv(sinew)sinew - kv(skold)skold] ∑ k)i-1 kv(sinew)sinew - kv(si-1new)si-1new
(39)
i+1
∑ Nkold[kv(skold)skold - kv(si-1new)si-1new]
Ni
new
)
k)i-1
kv(sinew)sinew - kv(si-1new)si-1new (40)
Birth/death terms can be calculated from population numbers. If further adaptation is needed after one mesh adaptation, the adaptation procedure will be repeated. Though the technique is difficult to implement, the obtained differential-algebraic equations and the strategy can be solved with standard integration software. Numerical Study To validate the proposed method, the seeded batch cooling crystallization data of potassium sulfate (K2SO4) is simulated. Jones and Mullin16 used the following assumptions (other assumptions can be found in their paper): well-mixed crystal suspension and no breakage and agglomeration of crystals. Also, the nucleation and growth kinetics found in their study are nonlinear functions of crystal size and process conditions as shown in the following.
B0 ) 2 × 108σ7.63
(41)
G ) 100σ2L0.5e-1/2400T
(42)
This study uses kinetic parameters found in their paper and the equation of the volumetric shape factor, kv, found in the paper of Budz et al.17
6234
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001
Table 1. Kernel Constants Obtained by Optimization used data set
Kagg
Kbre × 10-5
SSQconstant
SSQnatural
SSQconstant + SSQnatural
constant and natural cooling natural cooling scenario constant cooling scenario
2.955 2.767 0.313
8.480 3.719 8.444
34.16 43.40 10.09
14.77 17.03 645.04
48.93 60.43 655.13
for 0 < L < 100 µm, kv ) 0.896 exp(0.168L0.5 - 8.234 × 10-3L) (43) for 100 < L < 2000 µm, kv ) 4.460 exp(-0.0797L0.5 + 6.76 × 10-4L) (44) Unfortunately, the above equations are discontinuous at L ) 100 µm. Because a discontinuity can cause serious instability in integration, the function is reformulated with a sigmoid function which does not affect the original functions but converts them into a continuous function.
kv ) (1 - θ)[eq 43] + θ[eq 44] θ)
exp[0.2(L - 100)] 1 + exp[0.2(L - 100)]
(45) (46)
Among the four cooling strategies of Jones and Mullin’s study,16 their natural and constant supersaturation cooling scenarios are simulated (Figure 3). The temperature profile is programmed with time for 180 min for both.16
Constant cooling: T(t) ) 70 - 45 exp[-0.008356(180 - t)] (47) Natural cooling:
T(t) ) 15 + 45 exp(-0.008356t) (48)
For comparison of the obtained results with the experimental ones, the average number density is calculated as
n j i,numerical )
Ni vi+1 - vi
(49)
and we estimate number densities for the sizes used in the experimental result. They are based on average number densities calculated by eq 49 and obtained by the E01RAF routine (NAG Fortran library) by computing an interpolating rational function. For computational speed, the Jacobian matrix that has a sparse matrix form is evaluated analytically. We used 10-10 and 10-5 as the absolute and relative tolerance for the integration, and, vref and Nref in eq 32 are given as the approximated maximum values of the crystal size and number. They are 2 × 10-3 m and 106 no./kg of H2O, respectively. Without Agglomeration and Breakage. For the comparison, a natural cooling case neglecting agglomeration and breakage was simulated for Jones and Mullin’s study.16 To compare the simulation method, the population balance equation of eq 4 is transformed into a set of population balance equations based on the firstorder backward difference method (BDM) as follows. This case used uniform meshes of 50 and 1000.
In the result from adaptive mesh method (AMM), the nucleated and seed crystals grow as the crystals move (Figure 4). Especially, the crystal from seed appears in only 1 mesh because the method does not show any numerical diffusion. The result means that the models without agglomeration and breakage cannot explain the experimental results. However, the result from BDM shows very bad numerical diffusion, and a 1000-mesh result shows less numerical diffusion than a 50-mesh result. In the case of 50 meshes, the numerical diffusion is so important that it compensates the error made on the model, and the 50-mesh result seems to give an almost similar figure with the experimental result. So, numerical diffusion can lead to results similar to those in the experiment even with a wrong model. Therefore, it is underlined that a wrong model associated with a wrong method can make the user misunderstand the obtained reliable result. With Agglomeration and Breakage. Agglomeration and breakage are introduced to the simulation model. When the mass of nuclei mesh is greater than 0.1% of the total mass, a new mesh is added before the nuclei mesh. To show that the proposed method is suitable for the simulation of the crystallization process including agglomeration and breakage, an optimization to identify the kernel constant of agglomeration and breakage from the experimental result is performed. Because it is very difficult to estimate the four parameters (nucleation, growth rate, agglomeration kernel, and breakage kernel) simultaneously from CSD, it is assumed that the two equations of nucleation and growth rate are correct, and the optimization is performed only for the two mechanisms of agglomeration and breakage. The integrator used, DISCo, is restarted with the first-order method after every adaptation. Nevertheless, the optimization did not consume so many hours because one simulation took 20-32 s on a personal computer with a P-III 450 MHz processor. The optimization method used was the Simplex method.18 The parameters are obtained based on a logarithmic sum-of-squares objective function, and the tolerance used is 0.01. P
SSQ )
(ln n(v)exp - ln n(v)sim)i2 ∑ i)1
(51)
Here, P is the number of data points and is 17 in this case. For each search iteration, an updated parameter set was supplied by the optimization routine, a simulation was performed during 180 min, and a new SSQ was computed. Table 1 shows the result obtained based on the data set for the natural and constant cooling scenarios. The CSD obtained is shown in Figure 5a. For the comparison, the results optimized from two data sets respectively are drawn in parts b and c of Figure 5 and are summarized in Table 1. Conclusions
d(Gn) G(v) n(v) - G(v-∆v) n(v-∆v) ) dv ∆v
(50)
A simulation method for the crystalization processes was generated including breakage, growth, and ag-
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001 6235
glomeration, and breakage is proposed. The method is based on the discretized equations combined with the characteristics method. It is a very stable and accurate method because the difficulties on numerical diffusion and stability are avoided. This study assumed binary breakage and introduced a formula to distribute the broken particles into smaller meshes conserving number and mass. The method of characteristics is combined with the adaptive mesh method in order to handle stiff nucleation problems. The combined method improves the accuracy and efficiency of numerical computations. The proposed adaptive mesh method uses the strategy to equidistribute the arc length mesh function. Also, this study suggested the strategy to conserve number and mass during mesh adaptation. An example for the potassium sulfate system showed that the method is suitable for the simulation of crystallization processes. Acknowledgment This work was supported by the Korea Science and Engineering Foundation (995-1100-001-2) and the French Foreign Office. Notation A ) agglomeration kernel B ) breakage kernel b ) breakage function Ba ) birth rate of particles due to agglomeration Bb ) birth rate of particles due to breakage B0 ) nucleation rate at the smallest measurable size c ) concentration Da ) death rate of particles due to agglomeration Db ) death rate of particles due to breakage F ) flow rate f ) mesh function G ) growth rate H ) enthalphy Kagg ) agglomeration kernel constant Kbre ) breakage kernel constant Kg ) growth rate constant Kn ) nucleation rate constant kv ) volume shape factor L ) particle size M ) total number of meshes n ) population number density of particles N ) total number of particles Q ) heat duty S ) nucleation rate of particles si ) representative size for ith mesh T ) temperature t ) time U ) molar holdup V ) suspension volume v ) volume of the particles v0 ) smallest measurable volume vi ) volume of the smallest particles in ith mesh vi+1 ) volume of the largest particles in ith mesh x ) liquid molar fraction y ) gas molar fraction
η ) defined by eq 23 Subscripts a ) agglomeration b ) breakage C ) crystallization g ) gas i ) component i j ) index k ) feed l ) liquid s ) solid
Literature Cited (1) Tavare, N. S. Industrial Crystallization; Plenum Press: New York, 1995. (2) Tavare, N. S.; Garside, J.; Chivate, M. R. Analysis of batch crystallizers. Ind. Eng. Chem. Process Des. Dev. 1985, 19, 653665. (3) Ramkrishna, D. The status of population balances. Rev. Chem. Eng. 1985, 3, 49-95. (4) Hounslow, M. J.; Ryall, R. L.; Marshall, V. R. A discretized population balance for nucleation, growth and aggregation. AIChE J. 1988, 34, 1821-1832. (5) Marchal, P.; Davis, R.; Klein, J. P.; Villermaux, J. Crystallization and precipitation engineeringsI. An efficient method for solving population balances in crystallization with agglomeration. Chem. Eng. Sci. 1990, 43, 59-67. (6) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretizationsIII. Nucleation, growth, and aggregation of particles. Chem. Eng. Sci. 1997, 52, 4659-4679. (7) Petzold, L. R. Observations on an adaptive moving grid method for one-dimensional systems of partial differential equations. Appl. Numer. Math. 1987, 3, 347-360. (8) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes; Academic Press: New York, 1971. (9) Hartel, R. W.; Randolph, A. D. Mechanisms and kinetic modeling of calcium oxalate crystal aggregation in a urinelike liquorsII. Kinetic modeling. AIChE J. 1986, 32, 1186-1195. (10) Zauner, R.; Jones, A. G. Determination of nucleation, growth, agglomeration and disruption kinetics from experimental precipitation data: the calcium oxalate system. Chem. Eng. Sci. 2000, 55, 4219-4232. (11) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretizationsI. A Fixed Pivot Technique. Chem. Eng. Sci. 1996, 51, 1311-1332. (12) Hill, P. J.; Ng, K. M. New discretization procedure for the breakage equation. AIChE J. 1995, 41, 1204-1216. (13) Li, S. Adaptive mesh methods and software for timedependent partial differential equations. Doctoral dissertation, The University of Minnesota, Minneapolis, MN, 1998. (14) Huang, W.; Ren, Y.; Russell, R. D. Moving mesh methods based on moving mesh partial differential equations. J. Comput. Phys. 1994, 113, 279-290. (15) Sargousse, A.; LeLann, J. M.; Jourda, L.; Joulia, X. DISCo: un nouvel environnement de simulation oriente´ objet. 2nd Francophone Conference on Modeling and Simulation, Oct 6-8, 1999, Nancy, France, 1999; pp 61-66. (16) Jones, A.; Mullin, J. Programmed cooling crystallization of potassium sulfate solutions. Chem. Eng. Sci. 1974, 29, 105118. (17) Budz, J.; Jones, A. G.; Mullin, J. W. On the shape-size dependence of potassium sulfate crystals. Ind. Eng. Chem. Res. 1987, 26, 820-824. (18) Edgar, T. F.; Himmelblau, D. M. Optimization of Chemical Processes; McGraw-Hill: New York, 1987.
Greek Letters δi,j ) Dirac delta function σ ) relative supersaturation F ) molar density ∆c ) supersaturation
Received for review May 15, 2001 Revised manuscript received October 2, 2001 Accepted October 8, 2001 IE010443R
