- Email: [email protected]
Numerical solutions of population balance models in preferential crystallization
Chemical Engineering Science 63 (2008) 1342 – 1352 www.elsevier.com/locate/ces
Numerical solutions of population balance models in preferential crystallization S. Qamar a,∗ , A. Ashfaq a,b , I. Angelov c , M.P. Elsner c , G. Warnecke a , A. Seidel-Morgenstern b,c a Institute for Analysis and Numerics, Otto-von-Guericke University, PSF 4120, 39106 Magdeburg, Germany b Institute for Process Engineering, Otto-von-Guericke University, PSF 4120, 39106 Magdeburg, Germany c Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstrasse 1, 39106 Magdeburg, Germany
Available online 26 July 2007
Abstract This article focuses on the implementation of numerical schemes to solve model equations describing preferential crystallization for enantiomers. Two types of numerical methods are proposed for this purpose. The first method uses high resolution finite volume schemes, while the second method is the so-called method of characteristics (MOC). On the one hand, the finite volume schemes which were derived for general system in divergence form are computationally efficient, give desired accuracy on coarse grids, and are robust. On the other hand, the MOC offers a technique which is in general a powerful tool for solving growth processes, has capability to overcome numerical diffusion and dispersion, gives highly resolved solutions, as well as being computationally efficient. Several numerical test examples for a preferential crystallization model with and without fines dissolution under isothermal and non-isothermal conditions are considered. The comparison of the numerical schemes demonstrates clear advantages of the finite volume schemes and the MOC for the current model. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Population balance models; Preferential crystallization of enantiomers; Fines dissolution; High resolution schemes; Method of characteristics; Nucleation and growth rates
1. Introduction Separation of chiral molecules is an important question in industry as many organic and inorganic molecules are chiral. Illustration of enantiomers using amino acids as an example is given in Fig. 1. Usually, one of the enantiomers shows the desired properties for therapeutic activities or metabolism, whereas the other enantiomer may be inactive or may even cause some undesired effects. Complementing the most commonly used classical resolution via formation of diastereomers, direct crystallization methods have become increasingly important in recent years. An attractive process is based on enantioselective preferential crystallization. This concept is usually applied to the so-called conglomerates, which form physical mixture of crystals where each crystal is enantiomerically pure. In solution such systems tend to reach an equilibrium ∗ Corresponding author. Tel.: +49 391 6712027; fax: +49 391 6718073.
E-mail address: [email protected] (S. Qamar). 0009-2509/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2007.07.045
state in which the liquid phase will have racemic (50/50) composition and the solid phase will consist of a mixture of crystals of both enantiomers. However, before approaching this state, it is possible to preferentially produce crystals of just one of the enantiomers after seeding with the corresponding homochiral crystals. The process is essentially based on the different crystal surface areas of both enantiomers provided initially. The principle of preferential crystallization processes can be illustrated in a ternary phase diagram, see Fig. 2. A saturated solution at temperature Tcryst + T is cooled down to Tcryst , where the cooled liquid should be in the metastable zone, i.e. zone in which any particle-free supersaturated solution will stay clear for a finite time, that means no spontaneous (primary) nucleation takes place. In Fig. 2 point A represents an initial mixture of two enantiomers and a solvent (e.g. L-threonine, D-threonine and water). An enantiomeric excess (as depicted) is beneficial for the separation process, but not strictly necessary. At this point the crystallizer is seeded with crystals of enantiomer E1 . As the solution is supersaturated, the seeds will grow and
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
Fig. 1. Illustration of enantiomers using amino acids as an example.
Solvent Metastable zone
Equilibrium point
Seeding with E1 seeds
Tcryst Solubility curves
E
M
A Tcryst + ΔT
E1
Real trajectories after seeding with E1
E2
Fig. 2. Principle of preferential crystallization illustrated in a ternary phase diagram.
secondary nucleation will occur. If hypothetically there would be no primary nucleation of enantiomer E2 , the crystallization process would end at point M, which is the equilibrium point for E1 only. But experimental results show that after some time primary nucleation of E2 is induced and therefore, after sufficiently long time, the trajectory is attracted by the common equilibrium point for E1 andE2 , i.e. point E. Population balance models (PBMs), introduced in the late seventies (Ramkrishna, 1979), are the widely used simulation tool for different processes encountered in several scientific and engineering disciplines. They can be used to describe the time evolution of one or more property distributions of a population of particles. PBMs are used to study e.g. precipitation, polymerization, crystallization, particle size distribution (PSD) of crushed material and rain drops, dispersed phase distributions in multiphase flows, and so on. Typical phenomena occurring in these fields include growth, nucleation, aggregation, fragmentation, and condensation, among others. However, in the present work only growth and nucleation processes will be considered.
1343
Since population balance equations (PBEs) can only be solved analytically for simplified cases, numerical schemes are usually needed. Many studies are therefore focused on the development of accurate and efficient numerical solutions of the PBEs, see Ramkrishna (1985, 2000). High resolution schemes, which were originally developed for gas dynamics, were also applied for the numerical solution of PBEs, see Ma et al. (2002a, b), Gunawan et al. (2004). Qamar et al. (2007) have also used high resolution finite volume schemes in order to solve one and two-dimensional PBEs describing crystallization processes. These high resolution schemes and those considered in the present work have already been used for the numerical solution of hyperbolic systems which arise in astrophysical flows, gas dynamics, detonation waves and multi-phase fluid flows. The aim of these schemes is to obtain desired high order accuracy on coarse grids and resolve sharp discontinuities to avoid numerical diffusion and numerical dispersion which leads to unphysical oscillations. Since these schemes are developed for general purposes, they can be applied to hyperbolic problems in conservation law form without knowing the details of physical characteristics. Kumar and Ramkrishna (1997) have combined their discretization technique on a non-uniform grid for aggregation and breakage terms with a method of characteristics (MOC) for growth terms in order to solve PBEs with simultaneous nucleation, growth and aggregation processes. They found that, in contrast to other numerical techniques, the MOC avoids the numerical dissipation error caused by the growth term discretization. For handling the nucleation term, which is usually difficult to treat, the MOC was combined with a procedure of adding a cell of the nuclei size at a given time level. Furthermore, Lim et al. (2002) applied high resolution spatial discretization methods (WENO schemes) and the MOC for the dynamic simulations of batch crystallization including nucleation, growth, aggregation and breakage kinetics. In the present work, a new model has been derived for preferential crystallization of enantiomers with fines dissolution. The model is further elaborated by considering the isothermal and non-isothermal conditions. We have implemented the high resolution scheme of Koren (1993) and the one given by LeVeque (2002) for the solution of PBEs in preferential crystallization processes. Furthermore, we have also implemented the MOC (Kumar and Ramkrishna, 1997; Lim et al., 2002). These methods are used for the first time to model such processes. The schemes are discrete in space while continuous in time. The resulting ordinary differential equations can be solved by any standard ODE solver. In this work an adaptive RK45 method is used, which is an embedded Runge–Kutta method of order four and five. The paper is organized as follows. In Section 2, we start with a description of the model for preferential crystallization with fines dissolution. In Section 3, we briefly derive our proposed numerical schemes for the solution of the PBEs considered. Section 4 includes some numerical test problems for the current model with and without fines dissolution and considers as well isothermal and non-isothermal cases. In Section 5, we give final conclusions and remarks.
1344
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
2. Batch preferential crystallization model with fines dissolution In this section, we provide a mathematical model for the simulation of batch preferential crystallization with fines dissolution. Fines dissolution might lead to bigger crystals and more narrow distributions. This can fulfill special requirements regarding the product quality and simplifies downstream processes like filtration. A simplified dynamic model of an ideally mixed batch crystallizer for isothermal and non-isothermal conditions is supplied. Here we are considering the simplest case, in which a recycle pipe is attached to the crystallizer and we assume that fines dissolve completely at the end of the pipe. Fig. 3 shows the schematic diagram of the preferential crystallizer considered with the pipe for fines dissolution. The population balance for the solid phase is given as jn(k) (x, t) jn(k) (x, t) 1 = −G(k) (t) − h(x)n(k) (x, t), jt jx 1
(2.1)
where k ∈ {p, c} is the superscript to denote the variables describing the preferred and counter enantiomer respectively. In Eq. (2.1) n(k) (x, t) and G(k) (t) represent the number density and growth rate of the corresponding enantiomer, respectively. In Eq. (2.1) 1 is the residence time for the crystallizer. It is defined as V1 1 = , V˙
(2.2)
where V1 is the volume of the crystallizer and V˙ is the volumetric flow rate. Again in Eq. (2.1), h(x) is a death function, which describes the dissolution of particles below some critical size. It can be defined as a step (Heaviside) function. Fig. 4 presents such a function with xcrit = 2 × 10−4 m. The mass balance for the liquid phase in the crystallizer can be described by the following equation: dm(k) (t) (k) (k) ˙ out (t) =m ˙ in (t) − m dt ∞ (k) x 2 n(k) (x, t) dx. − 3kv G (t)
(2.3)
0
kv is the volumetric factor, is the crystal density and Here, ∞ 2 x n(x, t) dx is the second moment of the PSD. 0 Because of the fines dissolution this equation has two mass (k) fluxes. The first one m ˙ out (t) is the mass flux in the liquid phase which is being taken out from the vessel due to outgoing (k) liquid. The second one m ˙ in (t) is the incoming mass flux in the crystallizer from the dissolution pipe. They are defined as (k) m ˙ out (t) = w (k) liq (T )V˙ (k)
(2.4)
(k)
m ˙ in (t) = m ˙ out (t − 2 ) kv + 1
∞
x 3 h(x)n(k) (x, t − 2 ) dx.
0
preferred crystals (seeded)
counter crystals (unseeded) pipe for fines dissolution
fines dissolution unit
annular settling zone
tank A Fig. 3. Single-batch process setup with fines dissolution.
(2.5)
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
1345 (c)
preferred enantiomer and weq is equilibrium mass fraction for the counter enantiomer. The mass fractions are defined as
1
w(k) := 0.6
m(k) , + m(c) + mW
k ∈ {p, c},
(2.11)
where mW is the mass of solvent water. The nucleation rate kinetics are described as (Angelov et al., 2007; Elsner et al., 2005)
xcrit
h (x)
m(p)
(p)
(p)
(p)
B (p) (t) = kb (S (p) (t) − 1)b 3 (t), (c)
B (c) (t) = kb e−b 0
x, [m]
Here w (k) represent mass fractions and 2 is the residence time of the dissolution unit (the pipe) which is defined as 2 =
V˙
,
(2.6)
where r is the radius and l is the length of the pipe. The birth of new particles (nucleation) is incorporated in this model as a left boundary condition for the PBE at minimum crystal size x = x0 : n(k) (x0 , t) =
B (k) (t) , G(k) (t)
k ∈ {p, c}.
(2.7)
The assumption behind this substitution is that if B is the rate of appearance of nearly zero-sized particles and if the growth does not depend on the size, then dN dN dx B(t) = = = n(x0 , t)G(t). (2.8) · dt x→x0 ,t dx dt x→x0 ,t Here N is the number of crystals. The given model covers the case of isothermal operation, i.e. the temperature is considered to be constant during the whole batch process. A natural extension of the process idea is to vary the temperature during the batch process according to a specified and meaningful temperature profile. The growth rate kinetics are described as (Angelov et al., 2007; Elsner et al., 2005) (k)
G (t) = kg (S
(k)
g
(t) − 1) .
(2.9)
Here S (k) denotes the relative supersaturation for the preferred and counter enantiomers S (k) (t) =
w (k) (t) (k) weq (t)
,
k ∈ {p, c}.
.
(2.12)
Note that the above model reduces to the case no fines dissolution when the last term on the right-hand side of (2.1) and the first two terms on the right-hand side of (2.3) are zero. Then Eqs. (2.4) and (2.5) are not needed.
Fig. 4. Death function h(x).
r 2 l
(c) / ln(S (c) (t))2
(2.10)
In Eq. (2.10), w (p) (t) is the mass fraction of the dissolved preferred enantiomer, w(c) (t) is the mass fraction of the dissolved (p) counter enantiomer, weq is equilibrium mass fraction for the
3. Numerical schemes In this section we introduce two types of numerical methods for the solution of the above PBMs in crystallization. The first type of method includes high resolution finite volume schemes, while the second method is the so-called method of characteristics. For the derivation and explanation of the numerical schemes, it is convenient to consider a single population balance equation (PBE) for a batch crystallization process. The application of the schemes to the desired model given above is then straightforward. We are mainly interested in seeded batch crystallization processes and assume that the nucleation takes place at minimum crystal size. In this case the PBE has the form jn(x, t) j[G(x, t)n(x, t)] = − + B0 (t)(x − x0 ) + Q(x, t), jt jx (3.1) where B0 (t) is the nucleation rate of particles of minimum crystal size x0 and is a Dirac delta distribution. For details of this model see Hounslow (1990). Furthermore, Q(x, t) is any source term such as the last term on the right hand side of (2.1). Since the number density outside the computational domain is assumed to be zero, the above formulation is equivalent to considering the following PBE and defining the ratio of nucleation and growth terms as a left boundary condition jn(x, t) j[G(x, t)n(x, t)] =− + Q(x, t), jt jx n(x0 , t) =
B0 (t) . G(t)
(3.2) (3.3)
3.1. High resolution schemes Here, we give the formulation of high resolution finite volume schemes for the numerical solutions of above PBE. For a detailed discussion on the schemes the reader is referred to our article (Qamar et al., 2006).
1346
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
xi -
xi +
1 2
1 2
x xi-1
xi
xi+1
Fig. 5. Cell centered finite volume grid.
In order to apply a numerical scheme, the first step is to discretize the computational domain which is the crystal length in the current study. Let N be a large integer, and denoted by (xi−1/2 )i∈{1,...,N +1} a mesh of [x0 , xmax ], where x0 is the minimum and xmax is the maximum crystal length of interest. As shown in Fig. 5, for each i = 1, 2, . . . , N, x represents the cell width, the points xi refer to the cell centers, and the points xi±1/2 represent the cell boundaries. The integration of (3.2) over the cell i = [xi−1/2 , xi+1/2 ] yields the following cell centered semidiscrete finite-volume scheme: jn(x, t) Q(x, t) dx. dx = −((Gn)x + − (Gn)x − ) + i i jt i i (3.4) Let ni := ni (t) and Qi := Qi (t) denote the average value of the number density and source term in each cell i , i.e. 1 1 n(x, t) dx, Qi (t) = Q(x, t) dx. (3.5) ni (t) = x i x i Using Eq. (3.5), Eq. (3.4) implies 1 jni = − ((Gn)x + − (Gn)x − ) + Qi i i jt x
for i = 1, 2, . . . , N, (3.6)
by Koren (1993) and uses the value = 13 , while the second scheme was given by LeVeque (2002) with = −1. Note that the high resolution schemes need flux limitations in order to preserve monotonicity. This avoids the possible occurrence of wiggles and negative solution values (Koren, 1993; Qamar et al., 2006). HR-=1/3 scheme: According to this scheme the flux at the right boundary xi+1/2 of the cell i in (3.6) is (Koren, 1993) (Gn)x + = Gx + (ni + 21 (ri+ )(ni − ni−1 )), i
i
where the flux limiting function according to Koren (1993) is defined as (ri+ ) = max(0, min(2ri+ , min( 13 + 23 ri+ , 2))).
ri+ =
ni+1 − ni +
. ni − ni−1 +
B0 (t) . Gx −
(Gn)x + = Gx + ni ,
(Gn)x + = Gx + (ni + 21 ( + i )(ni+1 − ni )),
(Gn)x − = Gx − ni−1 . i
i
(3.7)
High order accuracy is easily obtained by piecewise polynomial interpolation. One can take for instance 1+ 1− (Gn)x + = Gx + ni + (ni+1 −ni ) + (ni − ni−1 ) , i i 4 4 (3.8) where ∈ [−1, 1]. Similarly one can write an expression for (Gn)x − . Here is a parameter that has to be chosen from the i indicated range. For = −1, one gets the second order accurate fully one-sided upwind scheme, and for = 1, the standard second order accurate central scheme. For all other values of ∈ [−1, 1], a weighted blend is obtained between the central scheme and the fully one-sided upwind scheme. In these studies we are interested in two types of flux limited high resolution schemes. The first one was introduced
(3.11)
This expression has to be evaluated with a small parameter, e.g.
= 10−10 , to avoid division by zero (Koren, 1993). Similarly, one can calculate the flux on the left boundary of the cell i . The nucleation is given as a left boundary condition n(x0 , t) =
i
(3.10)
Here the argument ri+ of this function is the so-called upwind ratio of two consecutive solution gradients
where N denotes the total number of mesh elements in the computational domain. The accuracy of finite volume discretizations is mainly determined by the way in which the cell-face fluxes are computed. Assuming that the flow is in positive x-direction, a first order accurate upwind scheme is obtained by taking the backward differences. i
(3.9)
(3.12)
1
HR- = −1 scheme: In this scheme the expression for the numerical flux can be obtained from Eq. (3.8) by taking = −1. The flux at the right boundary xi+1/2 of the cell i in (3.6) is (LeVeque, 2002) i
i
(3.13)
where + i =
ni − ni−1 . ni+1 − ni
Here the limiting function is given by the van Leer flux limiter (van Leer, 1985) (r) =
|r| + r 1 + |r|
(3.14)
with the same boundary condition as given by (3.12). Note that in the above schemes ni are the cell averaged values as given by (3.5), while the Gx + are the values of the size i dependent growth at the right face of each cell i. As discussed in the book of LeVeque (2002), this selection for the growth term will lead to a formally second order accurate scheme.
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
A disadvantage is that they cannot be applied straightforward to include boundaries. In order to avoid this problem one has to use the first order (3.7) in the first two cells on the left-hand side and the last cells on the right-hand side of the computational domain. In case of negative growth, i.e. G < 0 (i.e. dissolution), the discretization points (3.7)–(3.13) have to be mirrored at the considered volume boundaries xi+ for i = 1, 2, . . . , N. Finally, the resultant system of ordinary differential equations (ODEs) can be solved by any standard ODEs solver. 3.2. Method of characteristics (MOC) For linear advection equations a good alternative to the above schemes is the MOC. For a scalar linear conservation law, for example PBE (3.2) in the present case, there exist characteristic curves along which information propagates. The mesh is moved with the characteristic speed, whereby the linear advection is treated exactly. This drastically reduces numerical diffusion in the solution of the scheme in comparison to other schemes which use some discretization technique for approximating the advection term. For the MOC we use the same discretization as given in the previous subsection for the finite volume schemes as initial mesh at time t = 0 and consider a moving mesh along characteristics with mesh poi nts xi+1/2 (t) for i = 0, 1, 2, . . . , N. Let us substitute the growth rate G(t, x) by (3.15)
Then Eq. (3.2) leads to j jn(x, t) + jt jx
dx · n(x, t) = Q(x, t). dt
(3.16)
Integration over the control volume i (t)=[xi−1/2 (t), xi+1/2 (t)] gives i (t)
where xi (t) = xi+1/2 (t) − xi−1/2 (t). After using the above definitions, Eq. (3.18) implies d [(xi+1/2 (t) − xi−1/2 (t))ni ] = xi (t) Qi . dt
(3.20)
By using the product rule and (3.15), the left-hand side of (3.20) gives d [(xi+1/2 (t) − xi−1/2 (t))ni ] dt dxi+1/2 dxi−1/2 dni = xi (t) + − ni dt dt dt
dni = xi (t) + Gi+1/2 − Gi−1/2 ni . dt
(3.21)
After replacing the left-hand side of (3.20) with (3.21) and dividing the resultant equation by xi (t) one gets dni ni = −(Gi+1/2 − Gi−1/2 ) + Qi . dt xi (t)
(3.22)
In summary to use MOC, we have to solve the following set of equations: dni ni = −(Gi+1/2 − Gi−1/2 ) + Qi , dt xi (t) dxi+1/2 = Gi+1/2 , ∀i = 1, 2, . . . , N dt
(3.23) (3.24)
with initial data
dx := G(x, t). dt
1347
xi+1/2 (t) jn(x, t) dx dx+ · n(x, t) = Q(x, t) dx. jt dt i (t) xi−1/2 (t)
n(xi , 0) = n0 (xi ).
(3.25)
Note that in case of size independent growth the first term on the right-hand side of (3.23) is zero. In order to incorporate nucleation into the algorithm, a new cell of nuclei size is added at a given time level. The total number of mesh points can be kept constant by deleting the last cell at the same time level. Hence, all the variables such as ni (t) and xi (t) are initiated at these time levels and the time integrator restarts. In this case the boundary condition is again the same as given by (3.12). The above system of ordinary differential equations can be solved by any standard ODE solver.
(3.17) The Leibniz formula (Abramowitz and Stegun, 1972) on the left-hand side of (3.17) implies
d n(x, t) dx = Q(x, t) dx. dt i (t) i (t)
(3.18)
Let again ni := ni (t) and Qi := Qi (t) denote, respectively, the average values of the number density and source term in each cell i . As in Eq. (3.5), they are defined as 1 n(x, t) dx, xi (t) i 1 Qi := Q(x, t) dx, xi (t) i
ni :=
(3.19)
4. Numerical test problem In order to validate our numerical schemes for the considered model of preferential crystallization with fines dissolution, we consider the following numerical test problem. In this test problem we have taken the minimum crystal size x0 = 0. The idea behind choosing this problem is based on practical considerations, see Angelov et al. (2007) and Elsner et al. (2005) and references therein. The initial number density function of the seeds of preferred enantiomer is given as 1 1 ln(x) − 2 (p) 1 n (x, 0) = √ , (4.1) · · exp − 2 ·
2 Ia x Ia =
kV · s (p) 3 (0). Ms
(4.2)
1348
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
Table 1 Parameters for the test problem Description
Symbols
Value
Growth rate constant Growth rate exponent Activation energy Nucleation rate constant (seeded) Nucleation rate exponent (seeded) Activation energy Nucleation rate constant (unseeded) Nucleation rate exponent (unseeded) Universal gas constant Density of crystals Density of water Volume shape factor Initial mass of p-enantiomer Initial mass of c-enantiomer Mass of water Mass of seeds Residence time in crystallizer Volume of crystallizer Residencetime in pipe Radius of pipe Length of pipe Constant (seeded) Constant (seeded) Constant (seeded) Constant (seeded) Constant (seeded) Constant (seeded) Constant (unseeded) Constant (unseeded) Constant (unseeded) Constant (unseeded) Constant (unseeded) Constant (unseeded)
kg,0 g EA,g (p) kb,0 b(p) EA,b (c) kb,0 b(c) R
4.62 × 10 1.0 75.6 × 103 3.24 × 1025 4.0 78.7 × 103 3.84 × 106 6.0 × 10−2 8.314 1250 1000 0.0248 0.100224 0.100224 0.799552 2.5 × 10−3 60 1.0 × 10−3 10 1.0 × 10−2 5.3 × 10−1 0.056157 0.00143838 −3.41777 × 10−6 0.00624436 −0.0039185 4.4174 × 10−5 0.0574049 0.0013584 −2.3638 × 10−6 0.0071391 −0.00383673 4.2345 × 10−5
w
kv mp (0) mc (0) mW Ms
1
V1
2
r l A0 A1 A2 B0 B1 B2 C0 C1 C2 D0 D1 D2
Crystals of the counter enantiomer are initially not present, i.e. n(c) (x, 0) = 0.
(4.3)
Here we assume = 0.3947, = −6.8263, while Ms is the mass of the initial seeds. The maximum crystal size that was expected is xmax =0.005 m. The interval [0, xmax ] is subdivided into 500 grid points. The final time for the simulation was taken as 600 min. The kinetic parameters considered in this problem are given in Table 1. They are taken to describe the crystallization of the enantiomers of the amino acid threonine in water. The temperature trajectory used for the simulation is as follows: T (t)[Co ] = − 1.24074 × 10−7 t 3 + 4.50926 × 10−5 t 2 − 4.05556 × 10−3 t + 33.
(4.4)
Here the time is taken in minutes. While considering this tem(p) perature trajectory, the constants kg and kb will become functions of temperature as given below: kg = kg,0 · e−EA,g /R(T +273.15) , (k)
(k)
kb = kb,0 · e−EA,b /R(T +273.15) .
Unit 8
(4.5)
m/min – kJ/mol 1/m3 min – kJ/mol 1/min – J/mol K kg/m3 kg/m3 – Kg Kg Kg Kg min m3 min m m – – – – – – – – – – – –
(k)
Here kg,0 ,EA,g ,kb,0 ,EA,b , and R are constants given in Table 1. The equilibrium mass fraction of one enantiomer depends on the mass fraction of the other enantiomer and on temperature. They are given as (p)
weq (t) =
2
T j (t)(Aj + Bj w (c) (t)),
(4.6)
T j (t)(Cj + Dj w (p) (t)).
(4.7)
j =0 (c) weq (t) =
2
j =0
The constants in above two equations are also given in Table 1. Note that the same Eqs. (4.6) and (4.7) are also used in the isothermal case, where the temperature remains constant. Figs. 6 and 7 show the number density plots for the preferred enantiomer for the isothermal condition at 33 ◦ C and for the non-isothermal cases, temperature follows Eq. (4.4). In this figure the number densities for the cases with and without fines dissolution are given. Each figure compares the different numerical methods used to solve the PBM. In the isothermal case with fines dissolution, the peak of number density resulting from nucleation is smaller in comparison to the peak of number density without fines dissolution. It is due to the fact
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
x 108 Without Fines Dissolution
x 108
First Order HR κ=-1 HR-κ=1/3 MOC
10 8
With Fines Dissolution
7 particle density [#/m]
particle density [#/m]
12
6 4 2
1349
First Order HR κ=-1 HR-κ=1/3 MOC
6 5 4 3 2 1 0
0 0
1
2
3
4 x 10
crystal size [m]
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
5
x 10-3
crystal size [m]
-3
Fig. 6. Comparison of the number density of preferred enantiomer for different numerical methods and isothermal case at 600 min.
x 109
x 109
Without Fines Dissolution
With Fines Dissolution
8
First Order HR-κ=-1 HR-κ=1/3 MOC
5 4
particle density [#/m]
particle density [#/m]
6
3 2 1
data 1 data 2 data 3 MOC
7 6 5 4 3 2 1
0
0 0
1
2
3
crystal size [m]
4
5 x 10-3
0
1
2
3
crystal size [m]
4
5 x 10-3
Fig. 7. Comparison of number density of preferred enantiomer for different methods and non-isothermal case at 600 min.
that fines are taken out from the crystallizer, heated in a recycle pipe, thus dissolved and sent back to the crystallizer. This increases the supersaturation of the solution in the crystallizer which consequently reduces the secondary nucleation. Hence the number density resulting from nucleation is also reduced. The number density resulting from initial seed distribution will increase in the case of fines dissolution because the supersaturation of solution in the crystallizer increases which will allow the nuclei to grow more in comparison to the case without fines dissolution. In the non-isothermal case, the temperature is a timedependent polynomial. As temperature is reduced the nucleation rate increases which in turn increases the number density resulting from nucleation. Hence in the non-isothermal case, where the temperature is reducing with time, the number density increases in comparison to the isothermal case. The number of peaks in the crystal size distribution depends on the temperature profile. The number density in the non-isothermal case with fines dissolution is better because the crystal size is larger than in the other cases. Fig. 8 shows the mass fraction plots for the preferred (p-) and counter (c-) enantiomer with and without fines dissolution. In the isothermal case, mass fraction of the p-enantiomer
decreases sharply because we seeded p-enantiomer and it crystallizes out. While c-enantiomer mass fraction stays constant at the beginning, it then decreases later because of spontaneous primary nucleation. At some time both curves will join, which is the point where both enantiomers have similar equilibrium levels. The curve for each enantiomer with fines dissolution is higher in comparison to the case without fines dissolution which is due to the fact that solubility increases in case of fines dissolution. In the non-isothermal case, the mass fractions are completely controlled by the temperature profile. Fig. 9 shows the supersaturation plots for both enantiomers. Their behavior is entirely dependent on temperature which is clear from Eqs. (4.6) and (4.7). Figs. 10–12 show the growth rate, nucleation rate and third moment plots for the H R − = −1 scheme. They are the same for the other schemes. At the beginning the growth and nucleation rate of the p-enantiomer are reduced significantly because of sharp changes in mass fraction and supersaturation. In the non-isothermal case, the nucleation rate for the c-enantiomer with fines dissolution achieves a very high nucleation rate value because primary nucleation requires a high level of supersaturation. The third moment plots in the isothermal case for the p-enantiomer stays constant at the end of the process, because of
1350
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
mass fraction wk [kg/kg]
p-without fines c-without fines p-with fines c-with fines
0.1 0.098 0.096 0.094
Non-isothermal Case mass fraction wk [kg/kg]
Isothermal Case 0.102
0.092 0
0.1 0.098 0.096 0.094 0.092 0.09 0.088 0.086 0.084
100 200 300 400 500 600
p-without fines c-without fines p-with fines c-with fines
0
time [min]
100 200 300 400 500 600 time [min]
Fig. 8. Comparison of mass fractions in the liquid phase for the preferred (p-) and counter (c-) enantiomer.
Isothermal Case p-without fines c-without fines p-with fines c-with fines
0.08 0.06 0.04 0.02
Non-isothermal Case 0.12 supersaturation Sk-1 [-]
supersaturation Sk-1 [-]
0.1
p-without fines c-without fines p-with fines c-with fines
0.1 0.08 0.06 0.04 0.02
0 0
100 200 300 400 500 600 time [min]
0
100 200 300 400 500 600 time [min]
Fig. 9. Comparison of the supersaturation for the preferred (p-) and counter (c-) enantiomer.
Fig. 10. Comparison of the growth rates for the preferred (p-) and counter (c-) enantiomer.
no further change in mass fraction and supersaturation. For the non-isothermal case, the observed behavior results from mass fraction, supersaturation and temperature profile. Generally, for all plots in the non-isothermal case, the trends of the plots will vary with a different temperature profiles and initial crystal size distribution. Tables 2 and 3 show the percentage errors in mass preservation for the cases without and with fines dissolution under isothermal and non-isothermal conditions, respectively. Both
tables show that in the case of finite volume schemes, there are negligible changes in the percentage errors when we increase the number of mesh points from N = 500 to 1000. However, with the MOC the percentage errors decreases as we increase the number of mesh points. The finite volume schemes have large percentage errors as compared to the MOC. In these tables we have also given the computational time for all schemes. It is clear from the tables that MOC requires less computational time as compared to the finite volume schemes.
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
x 10-4
Isothermal Case p-without fines c-without fines p-with fines c-with fines
4000
nucleation rate [m/min]
nucleation rate [m/min]
5000
3000 2000 1000 0 0
4 3.5 3 2.5 2 1.5 1 0.5 0 -0.5
100 200 300 400 500 600
1351
Non-isothermal Case p-without fines c-without fines p-with fines c-with fines
0
100 200 300 400 500 600
time [min]
time [min]
Fig. 11. Comparison of the nucleation rates for the preferred (p-) and counter (c-) enantiomer.
x 10-4
x 10-4
Isothermal Case 7 third moment μ3 [m3]
third moment μ3 [m3]
4 3 2 p-without fines c-without fines p-with fines c-with fines
1
Non-isothermal Case p-without fines c-without fines p-with fines c-with fines
6 5 4 3 2 1 0
0 0
100
200
300
400
500
600
0
100
200
time [min]
300
400
500
600
time [min]
Fig. 12. Comparison of third moments for the preferred (p-) and counter (c-) enantiomer. Table 2 Percentage errors in mass preservation for the schemes without fines dissolution Method
First order scheme HR- = −1 scheme HR- = 13 scheme MOC
Isothermal
Non-isothermal
CPU time (s) (isothermal)
N = 500
N = 1000
N = 500
N = 1000
N = 500
N = 1000
3.737 3.811 3.813 2.604
3.775 3.813 3.814 1.844
4.460 4.733 4.736 3.792
4.669 4.736 4.737 2.917
1.5 2.2 2.3 0.34
3.1 4.4 4.6 0.41
Table 3 Percentage errors in mass preservation for the schemes with fines dissolution Method
First order scheme HR- = −1 scheme HR- = 13 scheme MOC
Isothermal
Non-isothermal
CPU time (s) (isothermal)
N = 500
N = 1000
N = 500
N = 1000
N = 500
N = 1000
2.801 2.873 2.875 1.823
2.838 2.875 2.876 1.30
2.841 2.962 2.965 2.055
2.904 2.965 2.967 1.086
2.3 3.1 3.5 0.39
5.5 7.5 7.7 0.71
However, for all schemes the computational times are less than 8 sec with N = 1000 mesh points. In Tables 2 and 3 we have only presented the computational time for the isothermal case
and was found to be the same for the non-isothermal case. All the computations were performed on computer with 1.73 GHz processor and 2 GB RAM. The programs are written in
1352
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
C programming language and were compiled with gcc 3.3.5 compiler using Suse Linux version 9.3 operating system. From the discussion in this section it is clear that computational efficiency of the numerical schemes is very important. Solution time is certainly a very important variable in the context of real-time optimization and control calculations but not so much in the context of off-line calculations. Therefore the current schemes, particularly the introduction of the MOC to deal with the hyperbolic nature of the growth term, are needed when accurate PBM solutions are needed to compute realtime control action for particulate processes. For instance, an efficient schemes is needed in a predictive control algorithms that require repeated model solutions for the computation of the control action through the solution of finite-time optimal control problems. A further discussion on the development and application of predictive-based strategies for control of PSD can be found in the recent articles of Shi et al. (2006, 2005). Although we are not considering control problems in this work, the current schemes can be used in this direction as well. 5. Conclusions In this article two semidiscrete high resolution finite volume schemes and a MOC are proposed for solving PBMs in preferential crystallization. For handling the nucleation term in the MOC, a cell of nuclei size is added at a given time level. The resultant systems of ODEs are then solved by a standard ODEsolver. In this work an adaptive RK45 method is used for this purpose, which is an embedded Runge–Kutta method of order four and five. Both high resolution schemes and the MOC worked very well for the numerical test cases presented in this article. However, it was found that the method of characteristics resolves better the number density distribution and give less error in the mass balances as compared to the finite volume schemes. Also the computational time of the finite volume schemes are larger as compared to the MOC. Thus, the use of the MOC for the growth processes considered in this paper is quite successful. However, due to their generality and robustness, high resolution schemes are very successful for modeling crystallization processes in general. They can be especially important when PBEs are coupled with computational fluid dynamics (CFD). In this case the models become more complicated and the solutions have strong discontinuities. In such applications high resolution schemes will be good candidates as the application of the MOC alone may not be possible due its several limitations. For example in the case of CFD model, the flux terms can be non-linear functions of conservative variables which are in turn function of time and external coordinates. Hence, MOC will be not applicable to the CFD model. However, one may use MOC for the internal property variables and FVM for the external coordinates. Furthermore, MOC is already adaptive due its built in moving mesh procedure and highresolution schemes can be easily combined with an adaptive mesh refinement procedure which are more important in
PBE–CFD applications. The use of adaptive mesh refinement will not only reduce the numerical errors in the solutions but also the overall computational time. References Abramowitz, M., Stegun, I.A., 1972. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover, New York. p. 11. Angelov, I., Raisch, J., Elsner, M.P., Seidel-Morgenstern, A., 2007. Optimal operation of enantioseparation by batch-wise preferential crystallization Chemical Engineering Science, in press, doi:10.1016/j.ces.2007.07.023. Elsner, M.P., Fernández Menéndez, D., Alonso Muslera, E., SeidelMorgenstern, A., 2005. Experimental study and simplified mathematical description of preferential crystallization. Chirality 17, 183–195. Gunawan, R., Fusman, I., Braatz, R.D., 2004. High resolution algorithms for multidimensional population balance equations. A.I.Ch.E. Journal 50, 2738–2749. Hounslow, M.J., 1990. A discretized population balance for continuous systems at steady-state. A.I.Ch.E. Journal 36, 106–116. Koren, B., 1993. A robust upwind discretization method for advection, diffusion and source terms. In: Vreugdenhill, C.B., Koren, B., (Eds.), Numerical Methods for Advection–Diffusion Problems, Notes on Numerical Fluid Mechanics, vol. 45. Vieweg Verlag, Braunschweig, 1993, pp. 117–138 (Chapter 5). Kumar, S., Ramkrishna, D., 1997. On the solution of population balance equations by discretization—III. Nucleation, growth and aggregation of particles. Chemical Engineering Science 52, 4659–4679. van Leer, B., 1985. Upwind-difference methods for aerodynamic problems governed by the Euler equations, Lectures in Applied Mathematics, vol. 22 (Part 2). American Mathematical Society, Providence, RI, 1985, pp. 327–336. LeVeque, R.J., 2002. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, New York. Lim, Y.I., Lann, J-M.L., Meyer, X.M., Joulia, X., Lee, G., Yoon, E.S., 2002. On the solution of population balance equation (PBE) with accurate front tracking method in practical crystallization processes. Chemical Engineering Science 57, 3715–3732. Ma, A., Tafti, D.K., Braatz, R.D., 2002a. High resolution simulation of multidimensional crystal growth. Industrial and Engineering Chemistry Research 41, 6217–6223. Ma, A., Tafti, D.K., Braatz, R.D., 2002b. Optimal control and simulation of multidimensional crystallization processes. Computers & Chemical Engineering 26, 1103–1116. Qamar, S., Elsner, M.P., Angelov, I., Warnecke, G., Seidel-Morgenstern, A., 2006. A comparative study of high resolution schemes for solving population balances in crystallization. Computers & Chemical Engineering 30, 1119–1131. Qamar, S., Ashfaq, A., Elsner, M.P., Angelov, I., Warnecke, G., Seidel-Morgenstern, A., 2007. Adaptive high resolution schemes for multidimensional population balances in crystallization processes. Computers & Chemical Engineering 31, 1296–1311. Ramkrishna, D., 1979. Statistical models of cell populations. Advances in Biochemical Engineering 11, 1–47. Ramkrishna, D., 1985. The status of population balances. Reviews Chemical Engineering 3 (1), 49–95. Ramkrishna, D., 2000. Population Balances: Theory and Applications to Particulate Systems in Engineering. Academic Press, San Diego, CA. Shi, D., Mhaskar, P., El-Farra, N.H., Christofides, P.D., 2005. Predictive control of crystal size distribution in protein crystallization. Nanotechnology 16, S562–S574. Shi, D., El-Farra, N.H., Li, M., Mhaskar, P., Christofides, P.D., 2006. Predictive control of particle size distribution in particulate processes. Chemical Engineering Science 61, 268–281.
Numerical solutions of population balance models in preferential crystallization S. Qamar a,∗ , A. Ashfaq a,b , I. Angelov c , M.P. Elsner c , G. Warnecke a , A. Seidel-Morgenstern b,c a Institute for Analysis and Numerics, Otto-von-Guericke University, PSF 4120, 39106 Magdeburg, Germany b Institute for Process Engineering, Otto-von-Guericke University, PSF 4120, 39106 Magdeburg, Germany c Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstrasse 1, 39106 Magdeburg, Germany
Available online 26 July 2007
Abstract This article focuses on the implementation of numerical schemes to solve model equations describing preferential crystallization for enantiomers. Two types of numerical methods are proposed for this purpose. The first method uses high resolution finite volume schemes, while the second method is the so-called method of characteristics (MOC). On the one hand, the finite volume schemes which were derived for general system in divergence form are computationally efficient, give desired accuracy on coarse grids, and are robust. On the other hand, the MOC offers a technique which is in general a powerful tool for solving growth processes, has capability to overcome numerical diffusion and dispersion, gives highly resolved solutions, as well as being computationally efficient. Several numerical test examples for a preferential crystallization model with and without fines dissolution under isothermal and non-isothermal conditions are considered. The comparison of the numerical schemes demonstrates clear advantages of the finite volume schemes and the MOC for the current model. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Population balance models; Preferential crystallization of enantiomers; Fines dissolution; High resolution schemes; Method of characteristics; Nucleation and growth rates
1. Introduction Separation of chiral molecules is an important question in industry as many organic and inorganic molecules are chiral. Illustration of enantiomers using amino acids as an example is given in Fig. 1. Usually, one of the enantiomers shows the desired properties for therapeutic activities or metabolism, whereas the other enantiomer may be inactive or may even cause some undesired effects. Complementing the most commonly used classical resolution via formation of diastereomers, direct crystallization methods have become increasingly important in recent years. An attractive process is based on enantioselective preferential crystallization. This concept is usually applied to the so-called conglomerates, which form physical mixture of crystals where each crystal is enantiomerically pure. In solution such systems tend to reach an equilibrium ∗ Corresponding author. Tel.: +49 391 6712027; fax: +49 391 6718073.
E-mail address: [email protected] (S. Qamar). 0009-2509/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2007.07.045
state in which the liquid phase will have racemic (50/50) composition and the solid phase will consist of a mixture of crystals of both enantiomers. However, before approaching this state, it is possible to preferentially produce crystals of just one of the enantiomers after seeding with the corresponding homochiral crystals. The process is essentially based on the different crystal surface areas of both enantiomers provided initially. The principle of preferential crystallization processes can be illustrated in a ternary phase diagram, see Fig. 2. A saturated solution at temperature Tcryst + T is cooled down to Tcryst , where the cooled liquid should be in the metastable zone, i.e. zone in which any particle-free supersaturated solution will stay clear for a finite time, that means no spontaneous (primary) nucleation takes place. In Fig. 2 point A represents an initial mixture of two enantiomers and a solvent (e.g. L-threonine, D-threonine and water). An enantiomeric excess (as depicted) is beneficial for the separation process, but not strictly necessary. At this point the crystallizer is seeded with crystals of enantiomer E1 . As the solution is supersaturated, the seeds will grow and
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
Fig. 1. Illustration of enantiomers using amino acids as an example.
Solvent Metastable zone
Equilibrium point
Seeding with E1 seeds
Tcryst Solubility curves
E
M
A Tcryst + ΔT
E1
Real trajectories after seeding with E1
E2
Fig. 2. Principle of preferential crystallization illustrated in a ternary phase diagram.
secondary nucleation will occur. If hypothetically there would be no primary nucleation of enantiomer E2 , the crystallization process would end at point M, which is the equilibrium point for E1 only. But experimental results show that after some time primary nucleation of E2 is induced and therefore, after sufficiently long time, the trajectory is attracted by the common equilibrium point for E1 andE2 , i.e. point E. Population balance models (PBMs), introduced in the late seventies (Ramkrishna, 1979), are the widely used simulation tool for different processes encountered in several scientific and engineering disciplines. They can be used to describe the time evolution of one or more property distributions of a population of particles. PBMs are used to study e.g. precipitation, polymerization, crystallization, particle size distribution (PSD) of crushed material and rain drops, dispersed phase distributions in multiphase flows, and so on. Typical phenomena occurring in these fields include growth, nucleation, aggregation, fragmentation, and condensation, among others. However, in the present work only growth and nucleation processes will be considered.
1343
Since population balance equations (PBEs) can only be solved analytically for simplified cases, numerical schemes are usually needed. Many studies are therefore focused on the development of accurate and efficient numerical solutions of the PBEs, see Ramkrishna (1985, 2000). High resolution schemes, which were originally developed for gas dynamics, were also applied for the numerical solution of PBEs, see Ma et al. (2002a, b), Gunawan et al. (2004). Qamar et al. (2007) have also used high resolution finite volume schemes in order to solve one and two-dimensional PBEs describing crystallization processes. These high resolution schemes and those considered in the present work have already been used for the numerical solution of hyperbolic systems which arise in astrophysical flows, gas dynamics, detonation waves and multi-phase fluid flows. The aim of these schemes is to obtain desired high order accuracy on coarse grids and resolve sharp discontinuities to avoid numerical diffusion and numerical dispersion which leads to unphysical oscillations. Since these schemes are developed for general purposes, they can be applied to hyperbolic problems in conservation law form without knowing the details of physical characteristics. Kumar and Ramkrishna (1997) have combined their discretization technique on a non-uniform grid for aggregation and breakage terms with a method of characteristics (MOC) for growth terms in order to solve PBEs with simultaneous nucleation, growth and aggregation processes. They found that, in contrast to other numerical techniques, the MOC avoids the numerical dissipation error caused by the growth term discretization. For handling the nucleation term, which is usually difficult to treat, the MOC was combined with a procedure of adding a cell of the nuclei size at a given time level. Furthermore, Lim et al. (2002) applied high resolution spatial discretization methods (WENO schemes) and the MOC for the dynamic simulations of batch crystallization including nucleation, growth, aggregation and breakage kinetics. In the present work, a new model has been derived for preferential crystallization of enantiomers with fines dissolution. The model is further elaborated by considering the isothermal and non-isothermal conditions. We have implemented the high resolution scheme of Koren (1993) and the one given by LeVeque (2002) for the solution of PBEs in preferential crystallization processes. Furthermore, we have also implemented the MOC (Kumar and Ramkrishna, 1997; Lim et al., 2002). These methods are used for the first time to model such processes. The schemes are discrete in space while continuous in time. The resulting ordinary differential equations can be solved by any standard ODE solver. In this work an adaptive RK45 method is used, which is an embedded Runge–Kutta method of order four and five. The paper is organized as follows. In Section 2, we start with a description of the model for preferential crystallization with fines dissolution. In Section 3, we briefly derive our proposed numerical schemes for the solution of the PBEs considered. Section 4 includes some numerical test problems for the current model with and without fines dissolution and considers as well isothermal and non-isothermal cases. In Section 5, we give final conclusions and remarks.
1344
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
2. Batch preferential crystallization model with fines dissolution In this section, we provide a mathematical model for the simulation of batch preferential crystallization with fines dissolution. Fines dissolution might lead to bigger crystals and more narrow distributions. This can fulfill special requirements regarding the product quality and simplifies downstream processes like filtration. A simplified dynamic model of an ideally mixed batch crystallizer for isothermal and non-isothermal conditions is supplied. Here we are considering the simplest case, in which a recycle pipe is attached to the crystallizer and we assume that fines dissolve completely at the end of the pipe. Fig. 3 shows the schematic diagram of the preferential crystallizer considered with the pipe for fines dissolution. The population balance for the solid phase is given as jn(k) (x, t) jn(k) (x, t) 1 = −G(k) (t) − h(x)n(k) (x, t), jt jx 1
(2.1)
where k ∈ {p, c} is the superscript to denote the variables describing the preferred and counter enantiomer respectively. In Eq. (2.1) n(k) (x, t) and G(k) (t) represent the number density and growth rate of the corresponding enantiomer, respectively. In Eq. (2.1) 1 is the residence time for the crystallizer. It is defined as V1 1 = , V˙
(2.2)
where V1 is the volume of the crystallizer and V˙ is the volumetric flow rate. Again in Eq. (2.1), h(x) is a death function, which describes the dissolution of particles below some critical size. It can be defined as a step (Heaviside) function. Fig. 4 presents such a function with xcrit = 2 × 10−4 m. The mass balance for the liquid phase in the crystallizer can be described by the following equation: dm(k) (t) (k) (k) ˙ out (t) =m ˙ in (t) − m dt ∞ (k) x 2 n(k) (x, t) dx. − 3kv G (t)
(2.3)
0
kv is the volumetric factor, is the crystal density and Here, ∞ 2 x n(x, t) dx is the second moment of the PSD. 0 Because of the fines dissolution this equation has two mass (k) fluxes. The first one m ˙ out (t) is the mass flux in the liquid phase which is being taken out from the vessel due to outgoing (k) liquid. The second one m ˙ in (t) is the incoming mass flux in the crystallizer from the dissolution pipe. They are defined as (k) m ˙ out (t) = w (k) liq (T )V˙ (k)
(2.4)
(k)
m ˙ in (t) = m ˙ out (t − 2 ) kv + 1
∞
x 3 h(x)n(k) (x, t − 2 ) dx.
0
preferred crystals (seeded)
counter crystals (unseeded) pipe for fines dissolution
fines dissolution unit
annular settling zone
tank A Fig. 3. Single-batch process setup with fines dissolution.
(2.5)
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
1345 (c)
preferred enantiomer and weq is equilibrium mass fraction for the counter enantiomer. The mass fractions are defined as
1
w(k) := 0.6
m(k) , + m(c) + mW
k ∈ {p, c},
(2.11)
where mW is the mass of solvent water. The nucleation rate kinetics are described as (Angelov et al., 2007; Elsner et al., 2005)
xcrit
h (x)
m(p)
(p)
(p)
(p)
B (p) (t) = kb (S (p) (t) − 1)b 3 (t), (c)
B (c) (t) = kb e−b 0
x, [m]
Here w (k) represent mass fractions and 2 is the residence time of the dissolution unit (the pipe) which is defined as 2 =
V˙
,
(2.6)
where r is the radius and l is the length of the pipe. The birth of new particles (nucleation) is incorporated in this model as a left boundary condition for the PBE at minimum crystal size x = x0 : n(k) (x0 , t) =
B (k) (t) , G(k) (t)
k ∈ {p, c}.
(2.7)
The assumption behind this substitution is that if B is the rate of appearance of nearly zero-sized particles and if the growth does not depend on the size, then dN dN dx B(t) = = = n(x0 , t)G(t). (2.8) · dt x→x0 ,t dx dt x→x0 ,t Here N is the number of crystals. The given model covers the case of isothermal operation, i.e. the temperature is considered to be constant during the whole batch process. A natural extension of the process idea is to vary the temperature during the batch process according to a specified and meaningful temperature profile. The growth rate kinetics are described as (Angelov et al., 2007; Elsner et al., 2005) (k)
G (t) = kg (S
(k)
g
(t) − 1) .
(2.9)
Here S (k) denotes the relative supersaturation for the preferred and counter enantiomers S (k) (t) =
w (k) (t) (k) weq (t)
,
k ∈ {p, c}.
.
(2.12)
Note that the above model reduces to the case no fines dissolution when the last term on the right-hand side of (2.1) and the first two terms on the right-hand side of (2.3) are zero. Then Eqs. (2.4) and (2.5) are not needed.
Fig. 4. Death function h(x).
r 2 l
(c) / ln(S (c) (t))2
(2.10)
In Eq. (2.10), w (p) (t) is the mass fraction of the dissolved preferred enantiomer, w(c) (t) is the mass fraction of the dissolved (p) counter enantiomer, weq is equilibrium mass fraction for the
3. Numerical schemes In this section we introduce two types of numerical methods for the solution of the above PBMs in crystallization. The first type of method includes high resolution finite volume schemes, while the second method is the so-called method of characteristics. For the derivation and explanation of the numerical schemes, it is convenient to consider a single population balance equation (PBE) for a batch crystallization process. The application of the schemes to the desired model given above is then straightforward. We are mainly interested in seeded batch crystallization processes and assume that the nucleation takes place at minimum crystal size. In this case the PBE has the form jn(x, t) j[G(x, t)n(x, t)] = − + B0 (t)(x − x0 ) + Q(x, t), jt jx (3.1) where B0 (t) is the nucleation rate of particles of minimum crystal size x0 and is a Dirac delta distribution. For details of this model see Hounslow (1990). Furthermore, Q(x, t) is any source term such as the last term on the right hand side of (2.1). Since the number density outside the computational domain is assumed to be zero, the above formulation is equivalent to considering the following PBE and defining the ratio of nucleation and growth terms as a left boundary condition jn(x, t) j[G(x, t)n(x, t)] =− + Q(x, t), jt jx n(x0 , t) =
B0 (t) . G(t)
(3.2) (3.3)
3.1. High resolution schemes Here, we give the formulation of high resolution finite volume schemes for the numerical solutions of above PBE. For a detailed discussion on the schemes the reader is referred to our article (Qamar et al., 2006).
1346
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
xi -
xi +
1 2
1 2
x xi-1
xi
xi+1
Fig. 5. Cell centered finite volume grid.
In order to apply a numerical scheme, the first step is to discretize the computational domain which is the crystal length in the current study. Let N be a large integer, and denoted by (xi−1/2 )i∈{1,...,N +1} a mesh of [x0 , xmax ], where x0 is the minimum and xmax is the maximum crystal length of interest. As shown in Fig. 5, for each i = 1, 2, . . . , N, x represents the cell width, the points xi refer to the cell centers, and the points xi±1/2 represent the cell boundaries. The integration of (3.2) over the cell i = [xi−1/2 , xi+1/2 ] yields the following cell centered semidiscrete finite-volume scheme: jn(x, t) Q(x, t) dx. dx = −((Gn)x + − (Gn)x − ) + i i jt i i (3.4) Let ni := ni (t) and Qi := Qi (t) denote the average value of the number density and source term in each cell i , i.e. 1 1 n(x, t) dx, Qi (t) = Q(x, t) dx. (3.5) ni (t) = x i x i Using Eq. (3.5), Eq. (3.4) implies 1 jni = − ((Gn)x + − (Gn)x − ) + Qi i i jt x
for i = 1, 2, . . . , N, (3.6)
by Koren (1993) and uses the value = 13 , while the second scheme was given by LeVeque (2002) with = −1. Note that the high resolution schemes need flux limitations in order to preserve monotonicity. This avoids the possible occurrence of wiggles and negative solution values (Koren, 1993; Qamar et al., 2006). HR-=1/3 scheme: According to this scheme the flux at the right boundary xi+1/2 of the cell i in (3.6) is (Koren, 1993) (Gn)x + = Gx + (ni + 21 (ri+ )(ni − ni−1 )), i
i
where the flux limiting function according to Koren (1993) is defined as (ri+ ) = max(0, min(2ri+ , min( 13 + 23 ri+ , 2))).
ri+ =
ni+1 − ni +
. ni − ni−1 +
B0 (t) . Gx −
(Gn)x + = Gx + ni ,
(Gn)x + = Gx + (ni + 21 ( + i )(ni+1 − ni )),
(Gn)x − = Gx − ni−1 . i
i
(3.7)
High order accuracy is easily obtained by piecewise polynomial interpolation. One can take for instance 1+ 1− (Gn)x + = Gx + ni + (ni+1 −ni ) + (ni − ni−1 ) , i i 4 4 (3.8) where ∈ [−1, 1]. Similarly one can write an expression for (Gn)x − . Here is a parameter that has to be chosen from the i indicated range. For = −1, one gets the second order accurate fully one-sided upwind scheme, and for = 1, the standard second order accurate central scheme. For all other values of ∈ [−1, 1], a weighted blend is obtained between the central scheme and the fully one-sided upwind scheme. In these studies we are interested in two types of flux limited high resolution schemes. The first one was introduced
(3.11)
This expression has to be evaluated with a small parameter, e.g.
= 10−10 , to avoid division by zero (Koren, 1993). Similarly, one can calculate the flux on the left boundary of the cell i . The nucleation is given as a left boundary condition n(x0 , t) =
i
(3.10)
Here the argument ri+ of this function is the so-called upwind ratio of two consecutive solution gradients
where N denotes the total number of mesh elements in the computational domain. The accuracy of finite volume discretizations is mainly determined by the way in which the cell-face fluxes are computed. Assuming that the flow is in positive x-direction, a first order accurate upwind scheme is obtained by taking the backward differences. i
(3.9)
(3.12)
1
HR- = −1 scheme: In this scheme the expression for the numerical flux can be obtained from Eq. (3.8) by taking = −1. The flux at the right boundary xi+1/2 of the cell i in (3.6) is (LeVeque, 2002) i
i
(3.13)
where + i =
ni − ni−1 . ni+1 − ni
Here the limiting function is given by the van Leer flux limiter (van Leer, 1985) (r) =
|r| + r 1 + |r|
(3.14)
with the same boundary condition as given by (3.12). Note that in the above schemes ni are the cell averaged values as given by (3.5), while the Gx + are the values of the size i dependent growth at the right face of each cell i. As discussed in the book of LeVeque (2002), this selection for the growth term will lead to a formally second order accurate scheme.
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
A disadvantage is that they cannot be applied straightforward to include boundaries. In order to avoid this problem one has to use the first order (3.7) in the first two cells on the left-hand side and the last cells on the right-hand side of the computational domain. In case of negative growth, i.e. G < 0 (i.e. dissolution), the discretization points (3.7)–(3.13) have to be mirrored at the considered volume boundaries xi+ for i = 1, 2, . . . , N. Finally, the resultant system of ordinary differential equations (ODEs) can be solved by any standard ODEs solver. 3.2. Method of characteristics (MOC) For linear advection equations a good alternative to the above schemes is the MOC. For a scalar linear conservation law, for example PBE (3.2) in the present case, there exist characteristic curves along which information propagates. The mesh is moved with the characteristic speed, whereby the linear advection is treated exactly. This drastically reduces numerical diffusion in the solution of the scheme in comparison to other schemes which use some discretization technique for approximating the advection term. For the MOC we use the same discretization as given in the previous subsection for the finite volume schemes as initial mesh at time t = 0 and consider a moving mesh along characteristics with mesh poi nts xi+1/2 (t) for i = 0, 1, 2, . . . , N. Let us substitute the growth rate G(t, x) by (3.15)
Then Eq. (3.2) leads to j jn(x, t) + jt jx
dx · n(x, t) = Q(x, t). dt
(3.16)
Integration over the control volume i (t)=[xi−1/2 (t), xi+1/2 (t)] gives i (t)
where xi (t) = xi+1/2 (t) − xi−1/2 (t). After using the above definitions, Eq. (3.18) implies d [(xi+1/2 (t) − xi−1/2 (t))ni ] = xi (t) Qi . dt
(3.20)
By using the product rule and (3.15), the left-hand side of (3.20) gives d [(xi+1/2 (t) − xi−1/2 (t))ni ] dt dxi+1/2 dxi−1/2 dni = xi (t) + − ni dt dt dt
dni = xi (t) + Gi+1/2 − Gi−1/2 ni . dt
(3.21)
After replacing the left-hand side of (3.20) with (3.21) and dividing the resultant equation by xi (t) one gets dni ni = −(Gi+1/2 − Gi−1/2 ) + Qi . dt xi (t)
(3.22)
In summary to use MOC, we have to solve the following set of equations: dni ni = −(Gi+1/2 − Gi−1/2 ) + Qi , dt xi (t) dxi+1/2 = Gi+1/2 , ∀i = 1, 2, . . . , N dt
(3.23) (3.24)
with initial data
dx := G(x, t). dt
1347
xi+1/2 (t) jn(x, t) dx dx+ · n(x, t) = Q(x, t) dx. jt dt i (t) xi−1/2 (t)
n(xi , 0) = n0 (xi ).
(3.25)
Note that in case of size independent growth the first term on the right-hand side of (3.23) is zero. In order to incorporate nucleation into the algorithm, a new cell of nuclei size is added at a given time level. The total number of mesh points can be kept constant by deleting the last cell at the same time level. Hence, all the variables such as ni (t) and xi (t) are initiated at these time levels and the time integrator restarts. In this case the boundary condition is again the same as given by (3.12). The above system of ordinary differential equations can be solved by any standard ODE solver.
(3.17) The Leibniz formula (Abramowitz and Stegun, 1972) on the left-hand side of (3.17) implies
d n(x, t) dx = Q(x, t) dx. dt i (t) i (t)
(3.18)
Let again ni := ni (t) and Qi := Qi (t) denote, respectively, the average values of the number density and source term in each cell i . As in Eq. (3.5), they are defined as 1 n(x, t) dx, xi (t) i 1 Qi := Q(x, t) dx, xi (t) i
ni :=
(3.19)
4. Numerical test problem In order to validate our numerical schemes for the considered model of preferential crystallization with fines dissolution, we consider the following numerical test problem. In this test problem we have taken the minimum crystal size x0 = 0. The idea behind choosing this problem is based on practical considerations, see Angelov et al. (2007) and Elsner et al. (2005) and references therein. The initial number density function of the seeds of preferred enantiomer is given as 1 1 ln(x) − 2 (p) 1 n (x, 0) = √ , (4.1) · · exp − 2 ·
2 Ia x Ia =
kV · s (p) 3 (0). Ms
(4.2)
1348
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
Table 1 Parameters for the test problem Description
Symbols
Value
Growth rate constant Growth rate exponent Activation energy Nucleation rate constant (seeded) Nucleation rate exponent (seeded) Activation energy Nucleation rate constant (unseeded) Nucleation rate exponent (unseeded) Universal gas constant Density of crystals Density of water Volume shape factor Initial mass of p-enantiomer Initial mass of c-enantiomer Mass of water Mass of seeds Residence time in crystallizer Volume of crystallizer Residencetime in pipe Radius of pipe Length of pipe Constant (seeded) Constant (seeded) Constant (seeded) Constant (seeded) Constant (seeded) Constant (seeded) Constant (unseeded) Constant (unseeded) Constant (unseeded) Constant (unseeded) Constant (unseeded) Constant (unseeded)
kg,0 g EA,g (p) kb,0 b(p) EA,b (c) kb,0 b(c) R
4.62 × 10 1.0 75.6 × 103 3.24 × 1025 4.0 78.7 × 103 3.84 × 106 6.0 × 10−2 8.314 1250 1000 0.0248 0.100224 0.100224 0.799552 2.5 × 10−3 60 1.0 × 10−3 10 1.0 × 10−2 5.3 × 10−1 0.056157 0.00143838 −3.41777 × 10−6 0.00624436 −0.0039185 4.4174 × 10−5 0.0574049 0.0013584 −2.3638 × 10−6 0.0071391 −0.00383673 4.2345 × 10−5
w
kv mp (0) mc (0) mW Ms
1
V1
2
r l A0 A1 A2 B0 B1 B2 C0 C1 C2 D0 D1 D2
Crystals of the counter enantiomer are initially not present, i.e. n(c) (x, 0) = 0.
(4.3)
Here we assume = 0.3947, = −6.8263, while Ms is the mass of the initial seeds. The maximum crystal size that was expected is xmax =0.005 m. The interval [0, xmax ] is subdivided into 500 grid points. The final time for the simulation was taken as 600 min. The kinetic parameters considered in this problem are given in Table 1. They are taken to describe the crystallization of the enantiomers of the amino acid threonine in water. The temperature trajectory used for the simulation is as follows: T (t)[Co ] = − 1.24074 × 10−7 t 3 + 4.50926 × 10−5 t 2 − 4.05556 × 10−3 t + 33.
(4.4)
Here the time is taken in minutes. While considering this tem(p) perature trajectory, the constants kg and kb will become functions of temperature as given below: kg = kg,0 · e−EA,g /R(T +273.15) , (k)
(k)
kb = kb,0 · e−EA,b /R(T +273.15) .
Unit 8
(4.5)
m/min – kJ/mol 1/m3 min – kJ/mol 1/min – J/mol K kg/m3 kg/m3 – Kg Kg Kg Kg min m3 min m m – – – – – – – – – – – –
(k)
Here kg,0 ,EA,g ,kb,0 ,EA,b , and R are constants given in Table 1. The equilibrium mass fraction of one enantiomer depends on the mass fraction of the other enantiomer and on temperature. They are given as (p)
weq (t) =
2
T j (t)(Aj + Bj w (c) (t)),
(4.6)
T j (t)(Cj + Dj w (p) (t)).
(4.7)
j =0 (c) weq (t) =
2
j =0
The constants in above two equations are also given in Table 1. Note that the same Eqs. (4.6) and (4.7) are also used in the isothermal case, where the temperature remains constant. Figs. 6 and 7 show the number density plots for the preferred enantiomer for the isothermal condition at 33 ◦ C and for the non-isothermal cases, temperature follows Eq. (4.4). In this figure the number densities for the cases with and without fines dissolution are given. Each figure compares the different numerical methods used to solve the PBM. In the isothermal case with fines dissolution, the peak of number density resulting from nucleation is smaller in comparison to the peak of number density without fines dissolution. It is due to the fact
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
x 108 Without Fines Dissolution
x 108
First Order HR κ=-1 HR-κ=1/3 MOC
10 8
With Fines Dissolution
7 particle density [#/m]
particle density [#/m]
12
6 4 2
1349
First Order HR κ=-1 HR-κ=1/3 MOC
6 5 4 3 2 1 0
0 0
1
2
3
4 x 10
crystal size [m]
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
5
x 10-3
crystal size [m]
-3
Fig. 6. Comparison of the number density of preferred enantiomer for different numerical methods and isothermal case at 600 min.
x 109
x 109
Without Fines Dissolution
With Fines Dissolution
8
First Order HR-κ=-1 HR-κ=1/3 MOC
5 4
particle density [#/m]
particle density [#/m]
6
3 2 1
data 1 data 2 data 3 MOC
7 6 5 4 3 2 1
0
0 0
1
2
3
crystal size [m]
4
5 x 10-3
0
1
2
3
crystal size [m]
4
5 x 10-3
Fig. 7. Comparison of number density of preferred enantiomer for different methods and non-isothermal case at 600 min.
that fines are taken out from the crystallizer, heated in a recycle pipe, thus dissolved and sent back to the crystallizer. This increases the supersaturation of the solution in the crystallizer which consequently reduces the secondary nucleation. Hence the number density resulting from nucleation is also reduced. The number density resulting from initial seed distribution will increase in the case of fines dissolution because the supersaturation of solution in the crystallizer increases which will allow the nuclei to grow more in comparison to the case without fines dissolution. In the non-isothermal case, the temperature is a timedependent polynomial. As temperature is reduced the nucleation rate increases which in turn increases the number density resulting from nucleation. Hence in the non-isothermal case, where the temperature is reducing with time, the number density increases in comparison to the isothermal case. The number of peaks in the crystal size distribution depends on the temperature profile. The number density in the non-isothermal case with fines dissolution is better because the crystal size is larger than in the other cases. Fig. 8 shows the mass fraction plots for the preferred (p-) and counter (c-) enantiomer with and without fines dissolution. In the isothermal case, mass fraction of the p-enantiomer
decreases sharply because we seeded p-enantiomer and it crystallizes out. While c-enantiomer mass fraction stays constant at the beginning, it then decreases later because of spontaneous primary nucleation. At some time both curves will join, which is the point where both enantiomers have similar equilibrium levels. The curve for each enantiomer with fines dissolution is higher in comparison to the case without fines dissolution which is due to the fact that solubility increases in case of fines dissolution. In the non-isothermal case, the mass fractions are completely controlled by the temperature profile. Fig. 9 shows the supersaturation plots for both enantiomers. Their behavior is entirely dependent on temperature which is clear from Eqs. (4.6) and (4.7). Figs. 10–12 show the growth rate, nucleation rate and third moment plots for the H R − = −1 scheme. They are the same for the other schemes. At the beginning the growth and nucleation rate of the p-enantiomer are reduced significantly because of sharp changes in mass fraction and supersaturation. In the non-isothermal case, the nucleation rate for the c-enantiomer with fines dissolution achieves a very high nucleation rate value because primary nucleation requires a high level of supersaturation. The third moment plots in the isothermal case for the p-enantiomer stays constant at the end of the process, because of
1350
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
mass fraction wk [kg/kg]
p-without fines c-without fines p-with fines c-with fines
0.1 0.098 0.096 0.094
Non-isothermal Case mass fraction wk [kg/kg]
Isothermal Case 0.102
0.092 0
0.1 0.098 0.096 0.094 0.092 0.09 0.088 0.086 0.084
100 200 300 400 500 600
p-without fines c-without fines p-with fines c-with fines
0
time [min]
100 200 300 400 500 600 time [min]
Fig. 8. Comparison of mass fractions in the liquid phase for the preferred (p-) and counter (c-) enantiomer.
Isothermal Case p-without fines c-without fines p-with fines c-with fines
0.08 0.06 0.04 0.02
Non-isothermal Case 0.12 supersaturation Sk-1 [-]
supersaturation Sk-1 [-]
0.1
p-without fines c-without fines p-with fines c-with fines
0.1 0.08 0.06 0.04 0.02
0 0
100 200 300 400 500 600 time [min]
0
100 200 300 400 500 600 time [min]
Fig. 9. Comparison of the supersaturation for the preferred (p-) and counter (c-) enantiomer.
Fig. 10. Comparison of the growth rates for the preferred (p-) and counter (c-) enantiomer.
no further change in mass fraction and supersaturation. For the non-isothermal case, the observed behavior results from mass fraction, supersaturation and temperature profile. Generally, for all plots in the non-isothermal case, the trends of the plots will vary with a different temperature profiles and initial crystal size distribution. Tables 2 and 3 show the percentage errors in mass preservation for the cases without and with fines dissolution under isothermal and non-isothermal conditions, respectively. Both
tables show that in the case of finite volume schemes, there are negligible changes in the percentage errors when we increase the number of mesh points from N = 500 to 1000. However, with the MOC the percentage errors decreases as we increase the number of mesh points. The finite volume schemes have large percentage errors as compared to the MOC. In these tables we have also given the computational time for all schemes. It is clear from the tables that MOC requires less computational time as compared to the finite volume schemes.
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
x 10-4
Isothermal Case p-without fines c-without fines p-with fines c-with fines
4000
nucleation rate [m/min]
nucleation rate [m/min]
5000
3000 2000 1000 0 0
4 3.5 3 2.5 2 1.5 1 0.5 0 -0.5
100 200 300 400 500 600
1351
Non-isothermal Case p-without fines c-without fines p-with fines c-with fines
0
100 200 300 400 500 600
time [min]
time [min]
Fig. 11. Comparison of the nucleation rates for the preferred (p-) and counter (c-) enantiomer.
x 10-4
x 10-4
Isothermal Case 7 third moment μ3 [m3]
third moment μ3 [m3]
4 3 2 p-without fines c-without fines p-with fines c-with fines
1
Non-isothermal Case p-without fines c-without fines p-with fines c-with fines
6 5 4 3 2 1 0
0 0
100
200
300
400
500
600
0
100
200
time [min]
300
400
500
600
time [min]
Fig. 12. Comparison of third moments for the preferred (p-) and counter (c-) enantiomer. Table 2 Percentage errors in mass preservation for the schemes without fines dissolution Method
First order scheme HR- = −1 scheme HR- = 13 scheme MOC
Isothermal
Non-isothermal
CPU time (s) (isothermal)
N = 500
N = 1000
N = 500
N = 1000
N = 500
N = 1000
3.737 3.811 3.813 2.604
3.775 3.813 3.814 1.844
4.460 4.733 4.736 3.792
4.669 4.736 4.737 2.917
1.5 2.2 2.3 0.34
3.1 4.4 4.6 0.41
Table 3 Percentage errors in mass preservation for the schemes with fines dissolution Method
First order scheme HR- = −1 scheme HR- = 13 scheme MOC
Isothermal
Non-isothermal
CPU time (s) (isothermal)
N = 500
N = 1000
N = 500
N = 1000
N = 500
N = 1000
2.801 2.873 2.875 1.823
2.838 2.875 2.876 1.30
2.841 2.962 2.965 2.055
2.904 2.965 2.967 1.086
2.3 3.1 3.5 0.39
5.5 7.5 7.7 0.71
However, for all schemes the computational times are less than 8 sec with N = 1000 mesh points. In Tables 2 and 3 we have only presented the computational time for the isothermal case
and was found to be the same for the non-isothermal case. All the computations were performed on computer with 1.73 GHz processor and 2 GB RAM. The programs are written in
1352
S. Qamar et al. / Chemical Engineering Science 63 (2008) 1342 – 1352
C programming language and were compiled with gcc 3.3.5 compiler using Suse Linux version 9.3 operating system. From the discussion in this section it is clear that computational efficiency of the numerical schemes is very important. Solution time is certainly a very important variable in the context of real-time optimization and control calculations but not so much in the context of off-line calculations. Therefore the current schemes, particularly the introduction of the MOC to deal with the hyperbolic nature of the growth term, are needed when accurate PBM solutions are needed to compute realtime control action for particulate processes. For instance, an efficient schemes is needed in a predictive control algorithms that require repeated model solutions for the computation of the control action through the solution of finite-time optimal control problems. A further discussion on the development and application of predictive-based strategies for control of PSD can be found in the recent articles of Shi et al. (2006, 2005). Although we are not considering control problems in this work, the current schemes can be used in this direction as well. 5. Conclusions In this article two semidiscrete high resolution finite volume schemes and a MOC are proposed for solving PBMs in preferential crystallization. For handling the nucleation term in the MOC, a cell of nuclei size is added at a given time level. The resultant systems of ODEs are then solved by a standard ODEsolver. In this work an adaptive RK45 method is used for this purpose, which is an embedded Runge–Kutta method of order four and five. Both high resolution schemes and the MOC worked very well for the numerical test cases presented in this article. However, it was found that the method of characteristics resolves better the number density distribution and give less error in the mass balances as compared to the finite volume schemes. Also the computational time of the finite volume schemes are larger as compared to the MOC. Thus, the use of the MOC for the growth processes considered in this paper is quite successful. However, due to their generality and robustness, high resolution schemes are very successful for modeling crystallization processes in general. They can be especially important when PBEs are coupled with computational fluid dynamics (CFD). In this case the models become more complicated and the solutions have strong discontinuities. In such applications high resolution schemes will be good candidates as the application of the MOC alone may not be possible due its several limitations. For example in the case of CFD model, the flux terms can be non-linear functions of conservative variables which are in turn function of time and external coordinates. Hence, MOC will be not applicable to the CFD model. However, one may use MOC for the internal property variables and FVM for the external coordinates. Furthermore, MOC is already adaptive due its built in moving mesh procedure and highresolution schemes can be easily combined with an adaptive mesh refinement procedure which are more important in
PBE–CFD applications. The use of adaptive mesh refinement will not only reduce the numerical errors in the solutions but also the overall computational time. References Abramowitz, M., Stegun, I.A., 1972. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover, New York. p. 11. Angelov, I., Raisch, J., Elsner, M.P., Seidel-Morgenstern, A., 2007. Optimal operation of enantioseparation by batch-wise preferential crystallization Chemical Engineering Science, in press, doi:10.1016/j.ces.2007.07.023. Elsner, M.P., Fernández Menéndez, D., Alonso Muslera, E., SeidelMorgenstern, A., 2005. Experimental study and simplified mathematical description of preferential crystallization. Chirality 17, 183–195. Gunawan, R., Fusman, I., Braatz, R.D., 2004. High resolution algorithms for multidimensional population balance equations. A.I.Ch.E. Journal 50, 2738–2749. Hounslow, M.J., 1990. A discretized population balance for continuous systems at steady-state. A.I.Ch.E. Journal 36, 106–116. Koren, B., 1993. A robust upwind discretization method for advection, diffusion and source terms. In: Vreugdenhill, C.B., Koren, B., (Eds.), Numerical Methods for Advection–Diffusion Problems, Notes on Numerical Fluid Mechanics, vol. 45. Vieweg Verlag, Braunschweig, 1993, pp. 117–138 (Chapter 5). Kumar, S., Ramkrishna, D., 1997. On the solution of population balance equations by discretization—III. Nucleation, growth and aggregation of particles. Chemical Engineering Science 52, 4659–4679. van Leer, B., 1985. Upwind-difference methods for aerodynamic problems governed by the Euler equations, Lectures in Applied Mathematics, vol. 22 (Part 2). American Mathematical Society, Providence, RI, 1985, pp. 327–336. LeVeque, R.J., 2002. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, New York. Lim, Y.I., Lann, J-M.L., Meyer, X.M., Joulia, X., Lee, G., Yoon, E.S., 2002. On the solution of population balance equation (PBE) with accurate front tracking method in practical crystallization processes. Chemical Engineering Science 57, 3715–3732. Ma, A., Tafti, D.K., Braatz, R.D., 2002a. High resolution simulation of multidimensional crystal growth. Industrial and Engineering Chemistry Research 41, 6217–6223. Ma, A., Tafti, D.K., Braatz, R.D., 2002b. Optimal control and simulation of multidimensional crystallization processes. Computers & Chemical Engineering 26, 1103–1116. Qamar, S., Elsner, M.P., Angelov, I., Warnecke, G., Seidel-Morgenstern, A., 2006. A comparative study of high resolution schemes for solving population balances in crystallization. Computers & Chemical Engineering 30, 1119–1131. Qamar, S., Ashfaq, A., Elsner, M.P., Angelov, I., Warnecke, G., Seidel-Morgenstern, A., 2007. Adaptive high resolution schemes for multidimensional population balances in crystallization processes. Computers & Chemical Engineering 31, 1296–1311. Ramkrishna, D., 1979. Statistical models of cell populations. Advances in Biochemical Engineering 11, 1–47. Ramkrishna, D., 1985. The status of population balances. Reviews Chemical Engineering 3 (1), 49–95. Ramkrishna, D., 2000. Population Balances: Theory and Applications to Particulate Systems in Engineering. Academic Press, San Diego, CA. Shi, D., Mhaskar, P., El-Farra, N.H., Christofides, P.D., 2005. Predictive control of crystal size distribution in protein crystallization. Nanotechnology 16, S562–S574. Shi, D., El-Farra, N.H., Li, M., Mhaskar, P., Christofides, P.D., 2006. Predictive control of particle size distribution in particulate processes. Chemical Engineering Science 61, 268–281.