4940
Ind. Eng. Chem. Res. 2010, 49, 4940–4947
An Efficient Numerical Method for Solving a Model Describing Crystallization of Polymorphs Shamsul Qamar,*,†,‡ Saima Noor,† and Andreas Seidel-Morgenstern‡ Department of Mathematics, COMSATS Institute of Information Technology, Park Road Chak Shahzad, Islamabad, Pakistan, and Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstrasse 1, 39106 Magdeburg, Germany
Polymorphism, in which a chemical compound exhibits different crystal forms or structures, has significant influence on the processing and storage of some crystalline powders in pharmaceutical industry. Different crystal structures, the so-called polymorphs, have different physical and chemical properties, such as crystal morphology, solubility, and color. These properties can have profound effect on the performance of products. This fact has motivated several researchers in this field to analyze, simulate, and control the crystallization of polymorphs. In this article, an efficient and accurate numerical method is introduced for solving a model describing crystallization of polymorphs. The proposed method has two parts. In the first part, a coupled system of ordinary differential equations of moments and the solute concentration is numerically solved in the time domain of interest. The resulting values are used to get the discrete growth and nucleation rates in the same time domain. In the second part, these discrete values along with the initial crystal size distribution (CSD) are used to construct the final CSD. The method is applied to a model describing crystallization of polymorphs of L-glutamic acid. For a validation, the results of the proposed technique are compared with those from the finite volume schemes. The numerical results demonstrate the potential of our scheme for the simulation of current model with high efficiency and accuracy. 1. Introduction Many crystalline compounds exhibit polymorphism, in which a chemical compound has the ability to contain more than one form or crystal structure. Batch crystallization process, used for the production of fine chemicals, is a process of solid formation from a supersaturated liquid. In most often situations, polymorphism phenomena is also observed in the solid crystalline materials. Different polymorphs have large variations in their physical properties, such as their shape, solubility, melting point, color, and chemical reactivity. Therefore, polymorphs can highly effect the functionality and properties, such as bioavailability, morphology, and purity, of many kinds of materials. For that reason, the crystallization of these polymorphs has to be controlled for improving the product quality and for the minimization of production costs. Recently, a rapid growth of interest has been observed in polymorphic crystallization.1-6 The modeling of a crystallization process enables one to achieve the desired goals and to investigate the effect of different operating conditions. The observed data can be used for optimization and control of the crystal size distribution (CSD). Population balance models (PBMs) have been widely used for modeling crystallization processes since the mid-1960s.7,8 Currently, PBMs are widely used for describing and controlling several particulate processes including comminution, crystallization, granulation, flocculation, combustion, and polymerization.9 Generally, population balance equations (PBEs) are partial integro-differential equations and have been applied in various engineering fields.9-16 The analytical solution of PBMs are only possible in simple situations, therefore numerical solutions are required in practical problems. However, an accurate simulation of the given PBM * To whom correspondence should be addressed. E-mail: qamar@ mpi-magdeburg.mpg.de. † COMSATS Institute Islamabad. ‡ Max Planck Institute Magdeburg.
can be challenging if the growth (dissolution) process is dominant and the distribution function has very narrow peaks or has sharp discontinuities. In such situations, the first order schemes can produce a large amount of numerical dissipation which normally grows during the simulation. A numerical dissipation in the scheme can be reduced by using very refined mesh at the cost of large computational time. However, a grid of very small mesh size can generate truncation errors in the numerical results. As a result, the numerical solution may converge to a wrong and smeared solution.17,18 Because of these reasons, much effort has been invested to develop appropriate numerical schemes for producing physically realistic solutions. The numerical methods for solving PBMs can be broadly divided into six classes, such as the method of moments,7,15,19-22 the method of characteristics,23-25 the method of weighted residualororthogonalcollocation,26 theMonteCarlosimulation,27,28 the finite difference schemes/discrete population balances,13,29 and the high resolution finite volume methods.14,10,25,30 Moreover, there also exist different approaches which are not based on the assumption of an ideally mixed fluid.31-34 Recently, Qamar et al.35 have derived an efficient Laplace transform based numerical method for solving a one-dimensional batch crystallization model with size-independent growth rate. The Laplace transform converted the given PBE to a linear ordinary differential equation (ODE) which was solved analytically by assuming that growth and nucleation rates were already known from the moment system of ODEs. Afterward, an inverse Laplace transform was used to get an expression for the actual CSD coupled with an implicit function of time variable. In this article, we implement the above numerical scheme for solving a model describing crystallization of polymorphs of metastable R-form and stable β-form crystals of L-glutamic acid.4,5,36 Moreover, an alternative procedure is introduced for the derivation of the proposed numerical method. In this derivation, instead of the Laplace transform, the method of characteristics and Duhamel’s principle are employed. Once
10.1021/ie9018353 2010 American Chemical Society Published on Web 04/15/2010
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
again, the main ingredients include initial data for the CSDs of polymorphs and the solute concentration, a coupled system of ODEs of moments, and the solute concentration, as well as an expression for the CSDs obtained by utilizing the method of characteristics and Duhamel’s principle. The method works as follows: given the initial solute concentration and moments, one can solve the ODE system of moments and the solute concentration. The solution of this system gives us discrete values of moments and the solute concentration in the time domain of interest and can be used to get the discrete values of growth and nucleation rates. The growth and nucleation rates along with the initial CSD are sufficient to construct the final CSD. For a validation, the numerical results of our method are also compared with those from the high resolution finite volume scheme (FVS) of Koren.37 It was found that, this technique is efficient, accurate, and easy to implement as compared to the FVS and other numerical methods used for solving this model.4,5 The article is arranged as follows. In section 2, a onedimensional model for batch crystallization of L-glutamic polymorphs is presented. In section 3, the proposed numerical method, named construction technique, is derived. In section 4, numerical test problems are presented. Finally, section 5 gives conclusions and remarks. 2. Model for Batch Crystallization of Polymorphs A model for the crystallization of polymorphs of metastable R-form and stable β-form crystals of L-glutamic acid is considered.4,5,36 The breakage and agglomeration processes are neglected in this model. In the one-dimensional case, the size of crystals is defined by a characteristic length x. The balance law for solid crystals is given as8,38 ∂n(k) ∂n(k) ) B(k)δ(x - x0), + G(k) ∂t ∂x n(k)(0, x) ) n(k) 0 (x),
(t, x) ∈ R+2
n(k)(t, x0) ) 0 ) n(k)(t, ∞),
(1)
k ) {R, β} (2)
G(R) )
∂n(k) ∂n(k) ) 0, + G(k) ∂t ∂x n (0, x) ) (k)
n(k) 0 (x),
2 (t, x) ∈ R+ ,
n (t, x0) ) (k)
B(k) 0 G(k)
,
k ) {R, β}
(3)
n (t, ∞) ) 0 (4) (k)
)
∫
∞ i (k) xn 0
dx i ) 0, 1, 2, 3, · · ·
dm 3 × 103 )F(k)κV(k)G(k)µ(k) 2 dt Fsolv k){R,β}
(6)
m(0) ) m0
(7)
∑
Here, Fsolv represents the density of solvent [kg/m3], F(k) is the density of k-form crystals [kg/m3], κV(k) is the volumetric shape
(
(β)
(β) G(β) ) K(β) - 1)g exp g (S
(8) (9)
-η(β) S -1 (β)
)
(β) (β) (β) B(β) ) k(β) - 1)µ(R) - 1)µ(β) b,1(S 3 + kb,2(S 3
(10) (11)
Here, S(k) ) m/m(k) sat , and msat is a saturation concentration [g/kg] defined as (k) (k) (k) 2 msat ) a(k) 1 T + a2 T + a 3
(12)
where T is the solution temperature [°C]. The kinetic (R) (R) parameters k(R) correspond to the nucleation b , kg , and kd -3 -1 rate [m s ], the growth rate [m/s], and the dissolution [m/s] rate of R-form crystals, respectively. Similarly, k(β) b, j and -3 -1 k(β) s ] and the g correspond to the jth nucleation rate [m growth rate [m/s] of β-form crystals, respectively. Moreover, η(β) is a dimensionless quantity, and g(k) is the growth exponential constant of k-form crystals. In the case of dissolution, that is, S(R) < 1, the growth rate (R) G and the nucleation rate B(R) in eqs 8 and 9 are negative. To account for variation in the crystal growth due to change in temperature, the following Arrhenius equation is used4,36
(
(k) k(k) g ) kg,0 exp -
)
E(k) g , R(T + 273)
k ∈ {R, β}
(13)
(k) where k(k) g, 0 represents the pre-exponential factor [m/s], Eg is the activation energy [J/mol] for the growth rate of k-form crystals, and R denotes the gas constant [J/(mol K)], respectively.
3. Derivation of the One-Dimensional Construction Technique Using the PBE (eq 1), the moments definition (eq 5), and equation for the solute concentration (eq 6), we get the 8-moment system of ODEs for both R- and β-form crystals coupled with an ODE for the solute concentration, dµ(k) 0 ) B(k) dt
(5)
A balance law for liquid phase yields an ODE of the solute concentration [g/kg] which is a sum of the concentrations of both polymorphs, that is, m ) m(R) + m(β), and is given as4,36
{
(R) k(R) - 1)g(R) if S(R) g 1 g (S (R) - 1) k(R) otherwise d (S
Analogously, for β-form crystals they are described as
The ith moment of the k-form CSD is defined as µ(k) i
(k) 3
(R) B(R) ) k(R) - 1)µ(R) g (S 3
(k)
where R+ :) (0,∞) and n represents the CSD of k-form crystals. Here, G(k) and B(k) represent the nucleation and growth rates of k-form crystals, and δ is the Dirac delta distribution. Note that eq 1 can also be written as
4941
factor (dimensionless) defined such that V ) κV x represents the total volume of k-form crystals [m3]. Moreover, 103 is a constant [g/kg] for ensuring consistent units. The power laws of crystal growth and nucleation rates for R-form crystals are given as4,36 (k)
dµ(k) i (k) ) iG(k)µi-1 + xi0B(k), dt
i ) 1, 2, 3,
(14)
k ) {R, β}
(15) dm 3 × 103 )F(k)κV(k)G(k)µ(k) 2 dt Fsolv k){R,β}
∑
(16)
The initial data for the above system are given as (k) µ(k) i (0) ) µi0 ,
m(0) ) m0,
k ) {R, β},
i ) 0, 1, 2, 3 (17)
4942
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
In the case of dissolution of metastable R-form crystals, that is, S(R) < 1, the nucleation rate B(R) in eq 9 is negative. Hence, according to eq 14, the number of crystals will be reduced whenever the CSD includes small-sized crystals (the case in practice), that is, the zeroth order moments µ(R) 0 will be reduced. In pharmaceutical practice, the objective in a crystallization in which multiple polymorphs occur is to have a final crystal population consisting of only one polymorph. Provided that a temperature profile can be realized keeping the solution supersaturated with respect to one polymorph and undersaturated with respect to the other polymorph, one polymorph will dissolve while the other one will grow. The above ODE-system can be solved by any standard ODEsolver. In this study we have used the RK45 method which is an embedded Runge-Kutta method of order four and five. The numerical solution of eqs 14-17 give us eight moments and the concentration in the given computational domain 0 e t e tmax. From these data, one gets the discrete values of the growth rate G(k) (by using eqs 8 and 10) and the nucleation rate B(k) (by using eqs 9 and 11) in the same time domain. The discrete values of growth and nucleation rates along with the initial CSDs for k-form crystals are sufficient to construct the final CSDs of both polymorphs. To get an expression for the CSD of k-form crystals, we solve the given PBEs in eq 1 by using the method of characteristics and Duhamel’s principle.39 This is a main step in the derivation of the proposed numerical method. Growth rates and nucleation rates are not explicit functions of time but are functions of supersaturation which varies with time. However, for the simplicity of notations, we will represent them as explicit functions of time whenever needed in the derivation of our numerical scheme as follows. Such representations are especially necessary in the integral terms obtained from Duhamel’s principle. To proceed with the solution of eq 1, we first solve the homogeneous part by applying the method of characteristics as follows (k)
(k)
∂n ∂n ) 0, + G(k) ∂t ∂x
n(k)(0, x) ) n(k) 0 (x)
(18)
(19)
Let dx/dt ) G(k), then the above equation can be rewritten as dn(k) ∂n(k) ∂n(k) + G(k) ) dt ∂t ∂x
(20)
After using eq 18 we finally obtain (k)
dn )0 dt
ξ(k)(t, x) ) x -
(21)
In other words k ) {R, β}
dτ(k) ) G(k), dt
(23)
(24)
k ) {R, β}
(25)
After solving these ODEs, we get τ(k) at the discrete points of the given time domain. Next, we solve the inhomogeneous PBE of the form ∂n(k) ∂n(k) + G(k) ) Q(k)(t, x), ∂t ∂x
n(k)(0, x) ) n(k) 0 (x)
After applying the Duhamel’s principle, we obtain
∫Q t
n(k) ) nH(k) +
(s, ξ(k)(t - s, x)) ds
(k)
0
(26)
Due to eq 1, Q(k) :) Q(k)(t, x) has the form Q(k) ) B(k)δ(x - x0)
(27)
By using eq 27 in eq 26, we obtain
∫B t
n(k)(t, x) ) nH(k)(t, x) +
(k)
0
(s) δ(ξ(k)(t - s, x) - x0) ds (28)
Let us define V(k): ) u(k)(s(k), t) )
∫
t
s(k)
G(k)(a) da
(29)
and ds(k) ∂ ) (k) (u(k))-1(V(k), t) (30) dV(k) ∂V
Therefore, we have ds(k) ∂ 1 ) ) (k) (u(k))-1(V(k), t) ) ∂ (k) (k) -1 (k) dV(k) ∂V u ((u ) (V , t), t) ∂s -1 (31) G(k)(s(k)) In light of the above derivations, eq 28 with (u(k))-1 :) (u(k))-1(V(k), t) and x˜ :) x - x0 becomes
∫
u(k)(t,t))0
u(k)(0,t)
B(k)((u(k))-1) δ(x˜ G(k)((u(k))-1) t G(k)(a) da) dV(k) (u(k))-1
∫
After using the definition of Dirac delta distribution, the above equation simplifies to n (t, x) ) (k)
k ) {R, β}
(s) ds ) x - τ(k)
(k)
τ(k)(0) ) 0,
(22)
where, the subscript ‘H’ is used for representing a homogeneous solution and c(k) is a constant of integration. After applying the initial conditions given in eq 17, we obtain (k) nH(k) ) n(k) 0 (ξ ),
t
0
where, τ(k) ) ∫t0G(k)(s) ds. As pointed out before, G(k) is represented as an explicit function of time inside the integral. To find τ(k), we have included the following two ODEs in the moment system given by eqs 14-17
n(k)(t, x) ) nH(k)(t, x) + nH(k) ) c(k),
∫G
s(k) ) (u(k))-1(V(k), t),
By using the chain rule, we obtain dn(k) dx ∂n(k) ∂n(k) + ) dt ∂t dt ∂x
where ξ(k) :) ξ(k)(t, x) is a location where the characteristic through (t, x) intersects the x- axis in the backward direction. It is given as
nH(k)(t, x)
{
B(k)(s(k)) (k) , + G(k)(s(k)) x˜ ∈ [0, u (0, t)] (32) 0, otherwise
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010 Table 1. Initial Seed Distribution Parameters for r-Form and β-Form Crystals symbols (R)
κ κ(β) σ(R)
F(0) > F(s ) > F(t) and F'(s ) < 0, (k)
values
units
symbols
values
units
2 × 10 2 × 1010 2
µm
jx jx(β) σ(β)
30 50 3
µm µm µm
10
(R)
Table 2. Model Parameters for Polymorphs of L-Glutamic Acid4,36 symbols
values
units
symbols
kb(R) (R) kg, 0 g(R) (R) Eg kd(R) (β) kb, 1 (β) kb, 2 (β) kg, 0 η(β) (β) g E(β) g
3.05 × 107 6.54 1.859 4.31 × 104 3.5 × 10-5 7.28 × 106 4.85 × 108 3.84 × 1022 0.78 1.047 1.76 × 105
m-3 s-1 m s-1 J mol-1 m s-1 m-3 s-1 m-3 s-1 m s-1 J mol-1
Fsolv F(R) ) F(β) κV(R) κV(β) a(R) 1 a(R) 2 a(R) 3 a(β) 1 a(β) 2 a(β) 3 b1 b2 b3
values 990 1540 0.48 0.031 8.437 3.032 4.564 7.644 1.165 0.622 37 8 0.1
× 10-3 × 10-2 × 10-3 × 10-1
units kg m-3 kg m-3 g kg-1 °C-2 g kg-1 °C-1 g kg-1 g kg-1 °C-2 g kg-1 °C-1 g kg-1 h h h
where, k ){R, β}. Note that, in the case of dissolution, G(k) is negative and hence u(k)(0, t) < 0. Therefore, the condition x˜ ∈ (0, u(k)(0, t)] in eq 32 implies a zero contribution of the nucleation term in the CSDs. The last step is to find s(k) :) (u(k))-1 in eq 32 needed for any x ∈ (0, u(k)(0, t)]. This can be obtained by finding a zero of the nonlinear function F(s(k)): )
∫
t
s(k)
G(k)(a) da - x˜,
(k)
0 0. If G(R) < 0 then u(R)(0, t) ) ∫t0G(R)(a) da < 0. Therefore, due to the condition x˜ ∈ [0, u(R)(0, t)] in eq 32, the nucleation term will not contribute in the final CSD of R-form crystals. In other words, this implies the lower case 0 in the second term of eq 32. Case 3. Finally, it is possible that R-form crystals initially grow slowly for a certain time and then start dissolving rapidly afterward. In this case, the growth rate is again a decreasing function of time having positive values initially and negative values afterward. As an example, see the growth rate of R-form crystals in Figure 5. In that figure a jump discontinuity can be observed when the positive growth rate abruptly changes to a negative one. Moreover, the integral function in eq 29 remains strictly monotone and becomes negative after starting a rapid dissolution of crystals. Thus, u(k)(0, t) < 0 implies the lower case in the second part of eq 32. The above three situations are commonly encountered in polymorphic crystallization processes. In numerical computa-
F'(s(k)): ) -G(k)(s(k)) (33)
Thus, for a given t and x ∈ (0, u(k)(0, t)] one can find s(k) by using Newton’s formula (k) sj+1 ) s(k) j -
F(s(k) j ) F'(s(k) j )
(34)
where j represents the iteration number in the Newton formula. Note that, the value of s(k) is different for both R- and β-form crystals. Since G(k) and B(k) are only available at discrete points of the time domain, a linear interpolation is needed for calculating G(k)(s(k)) and B(k)(s(k)) with s(k) ∈ [0, t]. Existence and Uniqueness of the Solution. In eq 30an inverse function has been introduced which is only possible if the original function is strictly monotone. In our case, the original function is an integral over the growth rate G(k) defined by (eq 29). The proof of monotonicity will guarantee the existence and uniqueness of a zero of the function in eq 33 and will ensure the general applicability of the technique. According to eq 8, the growth rate G(R) can take positive or negative signs, depending on S(R) > 1 or S(R) e 1. As a result, the integral function in eq 29 may not be strictly monotone. However, in the polymorphic crystallization process the signs of growth rate are not changing repeatedly. Therefore, three cases predominantly occur in polymorphic crystallization practice: Case 1. As one possibility, both R-form and β-form crystals may grow, that is, G(k) > 0. In that case, it is easy to show that F(s(k)) in eq 33 has a unique non-negative root s(k) ∈ (0, t)
∫G t
0
(k)
(a) da g x˜
(35)
In other words, one can find s(k) ∈ [0, t] such that F(s(k)) ) 0 in eq 33. Because of eq 35, eq 33 implies
Figure 1. Problem 1: CSDs of R-form and β-form crystals at t ) 4 h.
4944
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
Figure 2. Problem 1: Zoomed plots of CSDs of R-form and β-form crystals at t ) 4 h.
tions it was found that Newton’s routine needs four to six iterations for getting a convergent solution with a tolerance of ε ) 10-4. 4. Numerical Results Problem 1: Both r-Form and β-Form Crystals Grow. In this problem, the proposed numerical method is applied to the given model describing crystallization of polymorphs of L-glutamic acid. For validation, the numerical results of our scheme are compared with those from the first order finite volume scheme and the high resolution scheme of Koren.30,37 The initial seeds distribution for both R- and β-form crystals are given by Gaussian distributions4,36 n(k)(0, x) )
κ(k)
√2πσ(k)
(
exp -
(x - jx(k))2 2(σ(k))2
)
(37)
With the choice of parameters given in Table 1, the distribution function is narrowly sharp enough to challenge a numerical scheme. The remaining parameters needed in the simulation are given in Table 2. The temperature profile in the numerical simulation is chosen as T(t) ) b1 + b2e-t/b3 - t
(38)
where time t is taken in the unit of hours and b1, b2, and b3 are given in Table 2. Moreover, this temperature profile is arbitrarily chosen for the numerical purpose and is not taken from experiments. A very similar temperature profile used by Hermanto et al.4,36 is based on experiments. The plots
for R- and β-form CSDs are shown in Figure 1. In Figure 2, the zoomed parts of the distributions coming from nucleations and seeds are given. It is clear from the figures that the numerical results of our proposed scheme are very accurate compared to those obtained from the first order upwind scheme and the Koren scheme. Moreover, zoomed plots show that all three schemes approximate the nucleation parts of distributions almost similarly. However, clear differences are visible in the parts of distributions coming from seeds. The results of first order scheme are very diffusive, and even the Koren scheme is also not capable of achieving the peaks of the resulting distributions. However, our proposed numerical scheme completely preserves the narrow peaks of distribution and is totally free from numerical dissipation and oscillations. Figure 3 compares the normalized moments obtained from our scheme with those obtained from the Koren scheme. In these plots, symbols are used for the Koren scheme and lines are used for our scheme. Both schemes give overlapping results. Figure 4 gives the solute concentration and the temperature profile as defined by eq 38. Finally, Table 3 gives a comparison of errors in the solute concentration balance and the computational cost. Our scheme gives less error in the concentration balance compared to the finite volume schemes. Moreover, our numerical scheme is even faster than the first order upwind scheme. The Koren scheme was found more expensive for this test problem. The program was written in Matlab 6.5 software under a Linux operating system and was compiled on a computer with Intel(R) Core 2 Duo processor with 2 GHz speed and 3.83 GB (RAM)
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
4945
Figure 4. Problem 1: Concentration and temperature profiles. Figure 3. Problem 1: Normalized moments of R-form and β-form CSDs. Koren scheme (symbols); constructed scheme (lines).
memory. The numerical results and above discussion verify the robustness and accuracy of our proposed numerical scheme. Problem 2: Dissolution of r-Form Crystals and Growth of β-Form Crystals. This problem is intended to show the capability of our numerical method in handling dissolution of metastable R-form crystals and growth of stable β-form crystals. Here, R-form crystals grow initially at a slower rate for a certain time and then starts dissolving rapidly, while β-form crystals grow during the whole simulation time. Once again, the same initial CSDs given by eq 36 are considered with different values of jx(k). Here, we take jx(R) ) 50 µm and jx(β) ) 60 µm. Moreover, instead of eq 37 we take a constant temperature of T ) 43 °C. For this value the saturated (R) (R) concentrations are given as msat ) 21.468g/kg and msat ) 15.75g/kg. Figure 5 shows the numerical results. In this figure, one can see the dissolution of R-form crystals and the growth rate of β-form crystals. The growth rate profiles show that R-form crystals grow initially at a very slow rate for certain time (t < 0.8 h) but then start dissolving at a faster rate, while the growth of β-form crystals remains positive during the whole simulation time. Finally, the concentration decreases rapidly when both type crystals are growing (t < 0.8 h) and its decreasing rate reduces when R-form crystals starts dissolving. Once again, the numerical results of our scheme are superior as compared to the upwind finite volumes schemes.
Table 3. Errors in Concentration Balance description
absolute error
relative error
constructed Koren first-order
-14
-16
1.42 × 10 6.7 × 10-3 6.65 × 10-2
6.14 × 10 2.93 × 10-4 2.87 × 10-3
CPU time (s) 1.2 817.2 4.98
5. Conclusions An efficient and accurate numerical scheme has been derived for solving a model describing batch crystallization of polymorphs. The main ingredients of the scheme are the initial CSDs of polymorphs and the solute concentration, a coupled ODE system of moments and the solute concentration, as well as the expressions of CSDs derived by using the method of characteristics and Duhamel’s principle. The proposed numerical method is highly efficient and accurate. Moreover, the computer implementation of the scheme is also very easy. For validation, the numerical results of our scheme were compared with those obtained from the finite volume schemes. The proposed numerical method has a potential to resolve narrow peaks of the given CSD very well compared to the high resolution finite volume schemes. Furthermore, the method has the capability to resolve sharp discontinuities with their correct positions. In the numerical test problems, our scheme completely preserves the peaks of distributions and numerical dissipations are completely absent. However, numerical dissipations are visible in the finite volume schemes results. It should be pointed out that no grid refinement in the crystal size domain is needed for our numerical method because we
4946
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
Figure 5. Problem 2: CSDs, concentration, and growths profiles.
are not solving any partial differential equation numerically. Instead, the CSD is directly obtained from an analytic expression. The accuracy of the CDS is mainly dependent on the accurate computation of the moment system and on the roots of a nonlinear function of time variable. It was observed that even a lower number of mesh elements gives a numerical solution with the same accuracy. Moreover, our scheme preserves the solute concentration balance up to the order of machine accuracy. Acknowledgment Partial support by the Higher Education Commission (HEC) of Pakistan is gratefully acknowledged. Literature Cited (1) Blagden, N.; Davey, R. Polymorphs take shape. Chem. Brit. 1999, 35, 44–47. (2) Brittain, H. G. The impact of polymorphism on drug development: A regulatory viewpoint. Am. Pharm. ReV. 2000, 3, 67–70. (3) Fujiwara, M.; Nagy, Z. K.; Chew, J. W.; Braatz, R. D. First-principles and direct design approaches for the control of pharmaceutical crystallization. J. Process Control 2005, 15, 493–504. (4) Hermanto, M. W.; Kee, N. C.; Tan, R. B. H.; Chiu, M. -S.; Braatz, R. D. Robust bayesian estimation of kinetics for the polymorphic transformation of L-glutamic acid crystals. AIChE J. 2008, 54, 3248–3259. (5) Ono, T.; Kramer, H. J. M.; ter Horst, J. H.; Jansens, P. J. Process modeling of the polymorphic transformation of L-glutamic acid. Cryst. Growth Des. 2004, 4, 1161–1167.
(6) Scho¨ll, J.; Bonalumi, D.; Vicum, L.; Mazzotti, M.; Mller, M. In situ monitoring and modeling of the solvent-mediated polymorphic transformation of L-glutamic acid. Cryst. Growth Des. 2006, 6, 881–891. (7) Hulburt, H. M.; Katz, S. Some problems in particle technology. Chem. Eng. Sci. 1964, 19, 555–574. (8) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes, 2nd ed.; Academic Press: New York, 1988. (9) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Systems in Engineering; Academic Press: New York, 2000. (10) Gunawan, R.; Fusman, I.; Braatz, R. D. High resolution algorithms for multidimensional population balance equations. AIChE J. 2004, 50, 2738–2749. (11) Hounslow, M. J.; Ryall, R. L.; Marshall, V. R. A discretized population balance for nucleation, growth, and aggregation. AIChE J. 1988, 34, 1821–1832. (12) Hounslow, M. J.; Wynn, E. J. W. Short-cut models for particulate processes. Comput. Chem. Eng. 1993, 17, 505–516. (13) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization. I. A fixed Pivot technique. Chem. Eng. Sci. 1996, 51, 1311–1332. (14) Ma, D. L.; Tafti, D. K.; Braatz, R. D. High-resolution simulation of multidimensional crystal growth. Ind. Eng. Chem. Res. 2002, 41, 6217– 6223. (15) Marchisio, D. L.; Vigil, R. D.; Fox, R. O. Quadrature method of moments for aggregation-breakage processes. J. Colloid Interface Sci. 2003, 258, 322–334. (16) McCoy, B. J. A population balance framework for nucleation, growth, and aggregation. Chem. Eng. Sci. 2002, 57, 2279–2285. (17) LeVeque, R. J. Finite Volume Methods for Hyperbolic Problems; Cambridge University Press: New York, 2002. (18) Qamar, S.; Ashfaq, A.; Warnecke, G.; Angelov, I.; Elsner, M. P.; Seidel-Morgenstern, A. Adaptive high-resolution schemes for multidimensional population balances in crystallization processes. Comput. Chem. Eng. 2007, 31, 1296–1311.
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010 (19) Barrett, J. C.; Jheeta, J. S. Improving the accuracy of the moments method for solving the aerosol general dynamic equation. J. Aerosol Sci. 1996, 27, 1135–1142. (20) Madras, G.; McCoy, B. J. Reversible crystal growth-dissolution and aggregation breakage: Numerical and moment solutions for population balance equations. Powder Technol. 2004, 143-144, 297–307. (21) McGraw, R. Description of aerosol dynamics by the quadrature method of moments. Aerosol Sci. Technol. 1997, 27, 255–265. (22) Motz, S.; Mitrovic´, A.; Gilles, E. D. Comparison of numerical methods for the simulation of dispersed phase systems. Chem. Eng. Sci. 2002, 57, 4329–4344. (23) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization. III. Nucleation, growth and aggregation of particles. Chem. Eng. Sci. 1997, 52, 4659–4679. (24) Lim, Y. I.; Lann, J. -M. L.; Meyer, X. M.; Joulia, X.; Lee, G.; Yoon, E. S. On the solution of population balance equation (PBE) with accurate front tracking method in practical crystallization processes. Chem. Eng. Sci. 2002, 57, 3715–3732. (25) Qamar, S.; Warnecke, G. Numerical solution of population balance equations for nucleation growth and aggregation processes. Comput. Chem. Eng. 2007, 31, 1576–1589. (26) Rawlings, J. B.; Witkowski, W. R.; Eaton, J. W. Modelling and control of crystallizers. Powder Technology 1992, 69, 3–9. (27) Smith, M.; Matsoukas, T. Constant-number Monte Carlo simulation of population balances. Chem. Eng. Sci. 1998, 53, 1777–1786. (28) Tandon, P.; Rosner, D. E. Monte Carlo simulation of particle aggregation and simultaneous restructuring. J. Colloid Interface Sci. 1999, 213, 273–286. (29) Kumar, J.; Peglow, M.; Warnecke, G.; Heinrich, S.; Mo¨rl, L. A discretized model for tracer population balance equation: Improved accuracy and convergence. Comput. Chem. Eng. 2006, 30, 1278–1292. (30) Qamar, S.; Elsner, M. P.; Angelov, I.; Warnecke, G.; SeidelMorgenstern, A. A comparative study of high resolution schemes for solving population balances in crystallization. Comput. Chem. Eng. 2006, 30, 1119– 1131.
4947
(31) John, V.; Roland, M.; Mitkova, T.; Sundmacher, K.; Tobiska, L.; Voigt, A. Simulations of population balance systems with one internal coordinate using finite element methods. Chem. Eng. Sci. 2009, 64, 733– 741. (32) Kulikov, V.; Briesen, H.; Grosch, L.; von Wedel, L.; Yang, A.; Marqurdt, W. Modular dynamic simulation for integrated particulate processes by means of tool integration. Chem. Eng. Sci. 2005, 60, 2069– 2083. (33) Wang, L.; Marchisio, D.; Vigil, R.; Fox, R. CFD simulation of aggregation and breakage processes in laminar Taylor-Couette flow. J. Colloid Interface Sci. 2005, 282, 380–396. (34) Woo, X. Y.; Tan, R. B. H.; Chow, P. S.; Braatz, R. D. Simulation of mixing effects in antisolvent crystallization using a coupled CFD-PDFPBE approach. Cryst. Growth Des. 2007, 6, 1291–1303. (35) Qamar, S.; Warnecke, G.; Elsner, M. P.; Seidel-Morgenstern, A. A Laplace transformation-based technique for reconstructing crystal size distribution regarding size-independent growth. Chem. Eng. Sci. 2008, 63, 2233–2240. (36) Hermanto, M. W.; Braatz, R. D.; Chiu, M. -S. High-order simulation of polymorphic crystallization using weighted essentially nonoscillatory methods. AIChE J. 2008, 55, 122–131. (37) Koren, B. A robust upwind discretization method for advection, diffusion and source terms. In Numerical Methods for AdVection-Diffusion Problems; Vreugdenhil, C. B., Koren, B., Eds.; Notes on Numerical Fluid Mechanics, Vol. 45; Vieweg Verlag: Braunschweig, 1993; Vol. 11, pp 7138. (38) Rawlings, J. B.; Miller, S. M.; Witkowski, W. R. Model identification and control of solution crystallization processes: A review. Ind. Eng. Chem. Res. 1993, 32, 1275–1296. (39) Mcowen, R. C. Partial differential Equations: Methods and Applications, 2nd ed.; Prentice Hall: NJ, 2002.
ReceiVed for reView November 19, 2009 ReVised manuscript receiVed March 31, 2010 Accepted April 08, 2010 IE9018353
