- Email: [email protected]
Numerical approximations of a population balance model for coupled batch preferential crystallizers
Applied Numerical Mathematics 59 (2009) 739–753 www.elsevier.com/locate/apnum
Numerical approximations of a population balance model for coupled batch preferential crystallizers S. Qamar a,∗ , I. Angelov b , M.P. Elsner b , A. Ashfaq a , A. Seidel-Morgenstern b , G. Warnecke a a Institute for Analysis and Numerics, Otto-von-Guericke University, PSF 4120, 39106 Magdeburg, Germany b Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstrasse 1, 39106 Magdeburg, Germany
Available online 1 April 2008
Abstract This article is concerned with the numerical approximations of population balance equations for modeling coupled batch preferential crystallization processes. The current setup consists of two batch crystallizers interconnected with two fines dissolution pipes. The crystallization of both enantiomers is assumed to take place in separate crystallizers after seeding with their corresponding crystals. The withdrawn fines are assumed to be dissolved in the dissolution unit after heating. crystallizer, the crystallizer temperature before entering to the opposite crystallizer. Two types of numerical methods are proposed for the simulation of this process. The first method uses high resolution finite volume schemes, while the second method is the so-called method of characteristics. 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 method of characteristics 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. A numerical test problem with both isothermal and non-isothermal conditions is considered here. The numerical results show clear advantages of the proposed schemes. © 2008 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Population balance models; Preferential crystallization of enantiomers; Coupled crystallizers; Fines dissolution; High resolution schemes; Method of characteristics; Nucleation and growth
1. Introduction Separation of chiral molecules is an important industrial process as many, especially organic, molecules are chiral. Usually, one of the enantiomers shows desired properties for therapeutic activities or metabolism, whereas the other enantiomer may be inactive or may even cause some undesired effects. An attractive and energy efficient way for separation of chiral substances is the enantioselective preferential crystallization. This concept is usually applied to the so-called conglomerates, i.e. a physical mixture of crystals where each crystal is enantiomerically pure. Such systems tend to reach an equilibrium state in which the liquid phase will have racemic composition and the solid * Corresponding author. Tel: +49 391 6712027; fax: +49 391 6718073.
E-mail address: [email protected] (S. Qamar). 0168-9274/$30.00 © 2008 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2008.03.033
740
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Fig. 1. Principle of preferential crystallization illustrated in ternary phase diagram.
phase will consist of a mixture of crystals of both enantiomers. However, before approaching this state, it is possible to preferentially produce just one of the enantiomers after seeding with homochiral crystals. The principle of a preferential crystallization process can be illustrated in a ternary-phase diagram. Fig. 1 depicts the solvent and two pure enantiomers in the corners of the diagram. Starting from a saturated solution at temperature Tcryst +T , the solution becomes supersaturated if it is rapidly cooled down to temperature Tcryst within the metastable zone. In this zone the liquid will remain free of particles for a limited amount of time, as no spontaneous (primary) nucleation will occur. In Fig. 1 point A represents the initial mixture of two enantiomers and solvent (e.g.: L-threonine, D-threonine and water). An enantiomeric excess (as depicted) is beneficial for the separation process, but not strictly necessary. If at point A the crystallizer is seeded with seeds of both enantiomers E1 and E2 , the two crystal types will start to grow, and induce secondary nucleation. The process will eventually end at point E, which represents the equilibrium point of the two enantiomers for the temperature Tcryst . However, this is not an example of enantioselective preferential crystallization, as the two enantiomers are being produced simultaneously. On the other hand, if at point A the crystallizer is seeded with homochiral seeds of type E1 only, the liquid phase composition tends towards point M. However, the experiments show that after certain time, spontaneous crystallization of the unseeded enantiomer (in this case E2 ) is observed, and therefore the trajectory is attracted again by the common equilibrium point E, see [3,4]. Therefore the crystallization conditions have to be designed carefully and the process must be stopped before a significant nucleation of the unseeded (counter) enantiomer occurs – at point B for example. Detailed treatises of the process of crystallization by entrainment can be found in the literature [2,7,17]. In this article, the crystallization of the two enantiomers is assumed to take place in two separate crystallizers, coupled by their fines dissolution loops, see Fig. 2. The extracted solution is screened by filters and assumed to be free of larger crystals. Therefore, only small particles are withdrawn to the fines dissolution loop. In order to assure a crystal-free liquid exchange, the withdrawn liquid in the fines dissolution loop is heated, so that the liquid becomes undersaturated and the withdrawn small particles dissolve. Before entering into the opposite crystallizer, this liquid is assumed to be cooled again. This ensures that there will be no negative effect on the particles inside the opposite crystallizer due to warm liquid flux. The attrition or agglomeration processes are not considered in the current study. Furthermore, it is assumed that crystals of the one enantiomer do not change the birth and nucleation rates of the counter enantiomer. Population balance models (PBMs) are the widely used simulation tool for different processes encountered in several scientific and engineering disciplines, see [15]. They can be used to describe the time evolution of one or more property distributions of individuals population. PBMs are used to study e.g., precipitation, polymerization, crystallization, particle size distribution (PSD) of crushed material and rain drops, dispersed phase size 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.
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
741
Fig. 2. Coupled batch process setup.
Since population balance equations (PBEs) can be solved analytically only for simplified cases, numerical schemes are usually needed. Many studies are therefore focused on the development of accurate and efficient numerical solutions of the population balance equations, see Ramkrishna [16]. High resolution schemes, which were originally developed for gas dynamics, were also applied for the numerical solution of PBEs, see Ma et al. [12], Gunawan et al. [5]. Qamar et al. [13,14] have also used high resolution finite volume schemes in order to solve one and two-dimensional population balance equations 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 [9] 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. [11] applied high resolution spatial discretization methods (WENO schemes) and the method of characteristics 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 coupled batch preferential crystallization process. The model is further elaborated by considering the isothermal and non-isothermal conditions. We have implemented the high resolution scheme of Koren [8] and that by LeVeque [10] for solving PBEs which simulate preferential crystallization process in two coupled crystallizers. Furthermore, we have also implemented the method of characteristics [9,11]. 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 RK4 method is used, which is a Runge–Kutta method of order four. The paper is organized as follows. In Section 2, we start with the model for two coupled batch preferential crystallizers which are interconnected with fines dissolution units. In Section 3, we briefly derive our proposed numerical schemes for solving the current PBEs. Section 4 includes a numerical test problem for the current model for isothermal and non-isothermal cases. In Section 5, we give conclusions and remarks. 2. Coupled preferential crystallizers model with fines dissolution In this section, we provide a mathematical model for the simulation of two batch preferential crystallizers which are interconnected with fines dissolution pipes. For more details about preferential crystallization the reader is referred to [3,4]. A simplified dynamic model of ideally mixed batch crystallizers for isothermal and non-isothermal conditions is simulated. Here, we are considering the simplest case, in which inflow and outflow pipes are attached to the crystallizers and we assume that fines dissolves completely at the end of the pipe. Furthermore, we are assuming the same volume and residence time for both crystallizers. Similar assumptions hold for the two dissolution pipes as their
742
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Fig. 3. Death function h(x).
volume and residence time are identical. The function of pipes (heat exchanger) is to provide sufficient heat to the fines so that they can dissolve completely. Fig. 2 shows the schematic diagram of the coupled preferential crystallization process. In vessel A the L-threonine is the preferred enantiomer, D-threonine is the counter enantiomer and in vessel B vice versa. The population balance equations describing the balance laws for the solid phase are as follows: (k)
(k)
1 ∂nα (x, t) ∂nα (x, t) = −G(k) − hα (x)n(k) α (t) α (x, t) ∂t ∂x τ1,α with initial data
(2.1)
(k)
n(k) α (x, 0) = n0,α (x),
(2.2)
where k ∈ {p, c} and α ∈ {A, B}. The notation k stands for preferred or counter enantiomer and α for crystallizer A or B. In Eq. (2.1) τ1,α is the residence time in the corresponding crystallizer. It is defined as Vα , (2.3) τ1,α = V˙α where Vα is the volume of each crystallizer and V˙α is the corresponding volumetric flux rate. In Eq. (2.1) hα (x) is a death function describing the dissolution of particles below some critical size. It can be defined as a step (Heaviside) function. Fig. 3 represents one possible death function which we have considered in our numerical test problems with xcrit = 2 × 10−4 m. The mass balances for the liquid phase in the crystallizers can be described by the following equations (k)
dmA (t) (k) =m ˙ (k) ˙ (k) out,A (t) − 3ρkv GA (t) in,A (t) − m dt
∞
x 2 n(k) A (x, t) dx,
(2.4)
0 (k) dmB (t)
dt
(k)
(k)
(k)
∞
=m ˙ in,B (t) − m ˙ out,B (t) − 3ρkv GB (t)
(k)
x 2 nB (x, t) dx
(2.5)
0
with initial data (k)
(k)
mA (0) = m0,A ,
(k)
(k)
mB (0) = m0,B .
(2.6)
Here, ρ is the crystal density and kv is the volume shape factor defined such that the volume of a crystal with length (k) x is kv x 3 . Because of the fines dissolution Eqs. (2.4) and (2.5) have four mass fluxes. The first one m ˙ in,A (t) is the (k)
incoming mass flux in the liquid phase into the crystallizer A. The second one m ˙ out,A (t) is the outgoing mass flux from (k)
(k)
crystallizer A into the dissolution pipe. The flux m ˙ in,B (t) is the incoming mass flux into crystallizer B and m ˙ out,B (t) is the outgoing flux from crystallizer B. They are defined as follows (k) ˙ m ˙ (k) out,α (t) = wα ρliq Vα ,
k ∈ {p, c}, α ∈ {A, B}, ∞ kv ρ (k) (k) m ˙ in,α (t) = m ˙ out,α (t − τ2 ) + x 3 hα (x)n(k) α (x, t − τ2 ) dx, τ1,α 0
(2.7) (2.8)
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
743
(k)
where wα are the mass fractions and ρliq is the liquid density. Here τ2 is the residence time of the dissolution unit (the pipe) interconnecting both crystallizers. It is defined as τ2 =
πr 2 l , V˙
(2.9)
where r and l are the radius and length of each pipe, respectively. Due to the time delay in (2.8), special attention is needed for the numerical approximation of this equation. To simplify our problem, τ2 is chosen to be an integer multiple of the fixed time step t. At each time step the values of (k) the mass flow rates to the dissolution pipes mout,α (t) are saved in vectors which can be then used in (2.8) for t − τ2 0. (k) Similarly, The values of nα (x, t) are saved in matrices for using in (2.8) for t − τ2 0. For the case t − τ2 < 0 we (k) take min,α (t) = 0. However, one can also use an adaptive ODE-solver and some interpolation formula to calculate the terms in (2.8). It should be stressed that at the beginning of the process, i.e. t = 0, the two tanks including pipes are treated as one unit. For this reason the initial mass fractions for both enantiomers are everywhere the same. The birth of new particles (nucleation) is incorporated in this model as a left boundary condition for the PDE at minimum crystal size x = x0 (k)
n(k) α (x0 , t) =
Bα (t) (k)
(2.10)
.
Gα (t)
The assumption for this substitution is that if Bα is the rate of appearance of nearly zero-sized particles and the growth rate is independent of size, then dNα dNα dx (k) · = = nα (x0 , t)G(k) (2.11) Bα (t) = α (t). dt x→x0 dx dt x→x0 Here Nα is the number of crystals formed in crystallizer α ∈ {A, B}. The growth rate kinetics are described as: (k) g G(k) α (t) = kg Sα (t) − 1 .
(2.12) (k)
Here kg is the growth rate constant and g is the growth rate exponent. The function Sα denotes the relative supersaturation for the corresponding enantiomer defined as: (k)
Sα(k) (t) =
wα (t) (k)
(2.13)
.
weq,α (t)
(k)
(k)
The term wα (t) is mass fraction of the dissolved enantiomer, while weq,α (t) is equilibrium mass fraction for the corresponding enantiomer. The mass fractions are defined as: (k)
wα(k) =
mα (p)
(c)
mα + mα + mW
(2.14)
,
where mW is the mass of solvent (here water). The nucleation rate kinetics are described as: b(p) (p) (p) (p) μ3,α (t), Sα (t) − 1
Bα(p) (t) = kb
Bα(c) (t) = kb(c) e
−
b(c) (c) ln(Sα (t))2
,
(2.15)
where kb(k) are the nucleation rate constants and b(k) are the nucleation rate exponents for k ∈ {p, c}. Here, μ(k) 3,α (t) (k)
represent the third moments of the number densities nα (x, t). The given model covers the case of isothermal operation, i.e. the temperature is kept constant during the batch process. A natural extension of the process would be to vary the temperature during the batch process according to a specified temperature profile.
744
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Fig. 4. Cell centered finite volume grid.
3. Numerical schemes In this section we introduce two types of numerical methods for the solution of the above population balance models 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 ∂[G(x, t)n(x, t)] ∂n(x, t) =− + B0 (t) δ(x − x0 ) + Q(x, t), (3.1) ∂t ∂x 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 [6]. 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 homogeneous PBE and defining the ratio of nucleation and growth terms as a left boundary condition, i.e. ∂[G(x, t)n(x, t)] ∂n(x, t) =− + Q(x, t), ∂t ∂x
with n(x0 , t) =
B0 (t) . G(t)
(3.2)
3.1. High resolution schemes Here, we give the formulation of high resolution finite volume schemes for the numerical solutions of the above PBE. For a discussion of the schemes the reader is referred to our article [14]. 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 denote by (xi− 1 )i∈{1,...,N +1} a mesh of [x0 , xmax ], where x0 2 is the minimum and xmax is the maximum crystal length of interest. As shown in Fig. 4, for each i = 1, 2, . . . , N , x represents the cell width, the points xi refer to the cell centers, and the points xi± 1 represent the cell boundaries. 2 The integration of (3.2) over the cell Ωi = [xi− 1 , xi+ 1 ] yields the following cell centered semidiscrete finite2 2 volume scheme for Fi± 1 := (Gn)i± 1 2 2 ∂n(x, t) dx = −(Fi+ 1 − Fi− 1 ) + Q(x, t) dx. (3.3) 2 2 ∂t Ωi
Ωi
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.4) ni (t) = x x Ωi
Ωi
Using Eq. (3.4), Eq. (3.3) implies 1 dni =− (F 1 − Fi− 1 ) + Qi , for i = 1, 2, . . . , N (3.5) 2 dt x i+ 2 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.
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Fx
i+ 21
= Fi ,
Fx
i− 21
= Fi−1 .
745
(3.6)
High order accuracy is easily obtained by piecewise polynomial interpolation. One can take for instance 1+κ 1−κ (3.7) (Fi+1 − Fi ) + (Fi − Fi−1 ), 4 4 where κ ∈ [−1, 1]. Similarly one can write an expression for Fi− 1 . Here κ is a parameter that has to be chosen from 2 the 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 by Koren [8] and uses the value κ = 13 , while the second scheme was given by LeVeque [10] 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 [8,14]. Fi+ 1 = Fi + 2
HR-κ = 1/3 scheme: According to this scheme the flux at the right boundary xi+ 1 of the cell Ωi in (3.5) is [8] 2
1 Fi+ 1 = Fi + Φ(ri+ 1 )(Fi − Fi−1 ), 2 2 2 where the flux limiting function Φ according to Koren [8] is defined as 1 2 + r 1,2 . Φ(ri+ 1 ) = max 0, min 2ri+ 1 , min 2 2 3 3 i+ 2
(3.8)
(3.9)
Here the argument ri+ 1 of this function is the so-called upwind ratio of two consecutive flux gradients 2
ri+ 1 = 2
Fi+1 − Fi + ε . Fi − Fi−1 + ε
(3.10)
This expression has to be evaluated with a small parameter, e.g. ε = 10−10 , to avoid division by zero [8]. 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) =
B0 (t) . Gx −
(3.11)
1
HR-κ = −1 scheme: In this scheme the expression for the numerical flux can be obtained from Eq. (3.7) by taking κ = −1. The flux at the right boundary xi+ 1 of the cell Ωi in (3.5) is [10] 2
1 Fi+ 1 = Fi + φ(θi+ 1 )(Fi+1 − Fi ), 2 2 2
where θi+ 1 = 2
Fi − Fi−1 . Fi+1 − Fi
(3.12)
Here the limiting function φ is given by the van Leer flux limiter [18] φ(θi+ 1 ) = 2
|θi+ 1 | + θi+ 1 2
2
1 + |θi+ 1 |
(3.13)
2
with the same boundary condition as given by (3.11). 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 scheme (3.6) 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.6)–(3.12) have to be mirrored at the considered volume boundaries xi+ 1 for i = 1, 2, . . . , N . Finally, the resultant 2 system of ordinary differential equations can be solved by any standard ODEs solver. In the current study we have used the RK4 method which is a Runge–Kutta method of order four. For approximating the integral terms in the mass balance equations (2.4) and (2.5) we have used the composite Simpson’s rule.
746
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
3.2. Method of characteristics (MOC) For linear advection equations a good alternative to the above schemes is the method of characteristics. 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 points xi+1/2 (t) for i = 0, 1, 2, . . . , N . Let us substitute the growth rate G(t, x) by dx := G(x, t). dt Then Eq. (3.2) leads to ∂n(x, t) ∂ dx + · n(x, t) = Q(x, t). ∂t ∂x dt
(3.14)
(3.15)
Integration over the control volume Ωi (t) = [xi− 1 (t), xi+ 1 (t)] gives 2
Ωi (t)
2
x 1 (t) i+ dx ∂n(x, t) dx + · n(x, t) 2 = Q(x, t) dx. ∂t dt x 1 (t) i− 2
Ωi (t)
The Leibniz formula [1] on the left-hand side of (3.16) gives d n(x, t) dx = Q(x, t) dx. dt Ωi (t)
(3.16)
(3.17)
Ωi (t)
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.4), they are defined as 1 1 n(x, t) dx, Qi := Q(x, t) dx, (3.18) ni := xi (t) xi (t) Ωi
Ωi
where xi (t) = xi+ 1 (t) − xi− 1 (t). After using (3.18), Eq. (3.17) gives 2
2
d xi+ 1 (t) − xi− 1 (t) ni = xi (t) Qi . (3.19) 2 2 dt By using the product rule and (3.14), the left-hand side of (3.19) gives dx 1 dxi− 1 dni d i+ 2 2 xi+ 1 (t) − xi− 1 (t) ni = xi (t) + − ni 2 2 dt dt dt dt dni + (Gi+ 1 − Gi− 1 )ni . = xi (t) (3.20) 2 2 dt After replacing the left-hand side of (3.19) with (3.20) and dividing the resultant equation by xi (t) one gets dni ni = −(Gi+ 1 − Gi− 1 ) + Qi . 2 2 dt xi (t) In summary to use MOC, we have to solve the following set of equations dni ni = −(Gi+ 1 − Gi− 1 ) + Qi , 2 2 xi (t) dt dxi+ 1 2 = Gi+ 1 , ∀i = 1, 2, . . . , N 2 dt
(3.21)
(3.22) (3.23)
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
747
with initial data n(xi , 0) = n0 (xi ). Note that, in case of size independent growth the first term on the right-hand side of (3.22) 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 the case of stiff nucleation a small time step or comparatively refined mesh in the vicinity of the nucleation is desirable to avoid the loss in mass and unphysical oscillations, see [11]. For using different mesh speeds near the left boundary one can use the technique introduced in [9]. Usually, there is no effect on the MOC when the mesh becomes more refined, except an increase of the computational time. Here, the boundary condition is again the same as given by (3.11). The above system of ordinary differential equations can be solved by any standard ODE solver. In the current study we have used the RK4 method which is a Runge–Kutta method of order four. For approximating the integral terms in the mass balance equations (2.4) and (2.5) we have used the composite Simpson’s rule. 4. Numerical test problems In order to validate our numerical schemes for the current model for preferential crystallization, we consider the following numerical test problem. In this problem we have taken the minimum crystal size x0 = 0. The initial number density of the preferred enantiomer in each crystallizer is given as 1 ln(x) − μ 2 kV ρ (p) (p) , Iα = · exp − μ (0). (4.1) nα (x, 0) = √ √ Ms,α 3,α x 2π σ Iα 2σ The corresponding number density of the counter enantiomer in each crystallizer is zero, i.e. n(c) α (x, 0) = 0.
(4.2)
We assume σ = 0.3947 m, μ = −6.8263, and Ms,α are the masses of the initial seeds. The maximum crystal size that was expected is xmax = 0.005 m which is subdivided into 500 grid points. The final time for the simulation was taken as 600 minutes. The kinetic parameters considered in this problem are given in Table 1. They are taken to describe the crystallization of the enantiomers of amino acid threonine in water. Note that enantiomer E1 is preferred in crystallizer A and enantiomer E2 is preferred in crystallizer B. The temperature trajectory used for the simulation is as follows T (t)[C o ] = −1.24074 × 10−7 t 3 + 4.50926 × 10−5 t 2 − 0.00405556t + 33. Here the time is taken in minutes. For this temperature trajectory, the constants kg and temperature as given below EA,g
kg = kg,0 e− R(T +273.15) ,
EA,b
kb = kb,0 e− R(T +273.15) . (k)
(k)
(4.3) (p) kb
will become functions of
(4.4)
(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 the temperature. They are given as (p)
weq,α (t) =
2
T i (t) Ai + Bi w (c) (t) ,
(4.5)
T i (t) Ci + Di w (p) (t) .
(4.6)
i=0 (c) weq,α (t) =
2 i=0
The constants in the above two equations are given in Table 1. Note that the same Eqs. (4.5) and (4.6) are also used in the isothermal case (the temperature remains constant). Fig. 5 shows the number density plots in tank A for the preferred enantiomer with isothermal 33 ◦ C and for the non-isothermal condition, temperature follows Eq. (4.3). Note that, E1 is the preferred enantiomer in Tank A and E2 is the preferred enantiomer in Tank B. Due to the same initial conditions (mass of seeds, initial CSD) and due to the phenomenon that enantiomers have identical properties in arhiral environments, the plots for the enantiomer E1 in
748
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Table 1 Parameters for test problem 1 and 2 Description
Symbols
Value
Unit
Growth rate constant Growth rate exponent Activation energy
kg,0 g EA,g
4.62 × 108 1.0 75.6 × 103
m/s – kJ/mol
Nucleation rate constant (seeded)
kb,0
(p)
3.24 × 1025
1/m3 s
Nucleation rate exponent (seeded) Activation energy
b(p)
4.0 78.7 × 103
– kJ/mol
Nucleation rate constant (unseeded)
kb,0
(c)
3.84 × 106
1/s
Nucleation rate exponent (unseeded) Universal gas constant Density of crystals Density of water Volume shape factor Initial mass of p-enantiomer
b(c)
6.0 × 10−2 8.314 1250 1000 0.0248 0.100224
– J/mol K kg/m3 kg/m3 – kg
Initial mass of c-enantiomer Mass of water Mass of seeds Residence time in crystallizer Volume of crystallizer Residence time 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)
m0,α mW Ms τ1 V1 τ2 r l A0 A1 A2 B0 B1 B2 C0 C1 C2 D0 D1 D2
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
kg kg kg min m3 min m m – – – – – – – – – – – –
EA,b
R ρ ρw kv (p) m0,α (c)
Table 2 Mass preservation in the schemes Method
Percentage error isothermal case
Percentage error non-isothermal case
First order scheme HR-κ = −1 scheme HR-κ = 1/3 scheme MOC
2.945 3.020 3.022 1.781
2.148 2.301 2.302 1.996
tank A and for the enantiomer E2 in tank B are the same. In the isothermal case, the first large peak in the number density distribution is due to high nucleation rate of the preferred enantiomer at the beginning of the process. While the second small number density peak is due to the initial crystals growth. In the non-isothermal case, the temperature is a decreasing function of time. Hence the peak in the number density generated by nucleation has a higher value as compared to the isothermal case because of the temperature profile. The number of peaks in the crystal size distribution for the non-isothermal case entirely depends on the temperature profile used. The small zigzags in the left-hand side plot of MOC is due to stiff nucleation which was treated as a boundary condition. In MOC the mesh moves to the right at each time step and we have to add one mesh cell at the left boundary to include the boundary conditions. In the case
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
749
Fig. 5. Comparison of the number density for the preferred enantiomer for different operation modes at t = 600 minutes.
Fig. 6. Three dimensional plots of number density for the preferred enantiomer by using MOC.
of stiff nucleation a small time step or comparatively refined mesh in the vicinity of the nucleation is desired to avoid such oscillations, see [9,11]. There is no visible difference between results of the given schemes, so one cannot say exactly which scheme is better one. Table 2 shows that percentage errors for mass preservation are minimum for the MOC as compared to the finite volume schemes. The reason can be the current stiff nucleation which is considered as a boundary condition in the finite volume schemes. In the case of high resolution schemes this stiff nucleation may produce over predictions. This may lead to the bigger errors in the mass balances. Furthermore, there is also no significant difference in mass fractions of the MOC and the finite volume schemes, as well as supersaturation, growth rate, nucleation rates and the third moment. Therefore we have only presented here the plots of the MOC. Figs. 6 and 7 show the three dimensional plot of number density for the preferred enantiomer and the counter enantiomer. Fig. 8 shows the mass fraction plots for the preferred (p-) and counter (c-) enantiomers. In the isothermal case, the mass fraction of the p-enantiomer decreases sharply because we have seeded the p-enantiomer and it crystallizes out. While c-enantiomer mass fraction stay 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. 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.5) and (4.6). Figs. 10 and 11 show the growth rate and nucleation rate plots. At the beginning, growth rate and nucleation rate for the p-enantiomer are reduced significantly because of sharp changes in mass fraction and supersaturation. In the non-isothermal case, nucleation rate for counter enantiomer achieves a very high value because of the high supersaturation level.
750
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Fig. 7. Three dimensional plots of number density for the counter enantiomer using MOC.
Fig. 8. Comparison of the mass fractions for the preferred (p-) and counter (c-) enantiomer using MOC.
Fig. 9. Comparison of the supersaturation for the preferred (p-) and counter (c-) enantiomer using MOC.
Fig. 12 shows the trajectories of the third moments. These plots in the isothermal case for the p-enantiomer stays constant at the end of the process, because of no further change in mass fraction and supersaturation. In the nonisothermal case the third moment does not approach a steady state because of the decreasing temperature in the crystallizer.
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Fig. 10. Comparison of growth rate for the preferred (p-) and counter (c-) enantiomer using MOC.
Fig. 11. Comparison of nucleation rate for the preferred (p-) and counter (c-) enantiomer using MOC.
Fig. 12. Comparison of the third moment for the preferred (p-) and counter (c-) enantiomer using MOC.
751
752
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
5. Conclusions In this article two high resolution finite volume schemes and a method of characteristics were proposed for solving population balance models in two coupled preferential crystallizers. For the process of preferential crystallization there are two main advantages of considering two coupled crystallizers which are interconnected by two fines dissolution units. The first one is that one gets both enantiomers at the same time in separate crystallizers. Secondly, because of the fines dissolution, the amount of small particles reduces which further enhances the particle growth. For handling the nucleation term in the method of characteristics, a cell of nuclei size is added at a given time level. The resultant system of ODEs are then solved by a standard ODE-solver. In this work we have used the RK4 method which is a Runge–Kutta method of order four. Both high resolution schemes and the method of characteristics worked very well for the numerical test problems which are considered in this article. However, it was found that the method of characteristics resolves better the number density distribution and give smaller errors 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 method of characteristics. The use of MOC for the current growth processes looks quite successful. However, due to their generality and robustness to model perturbations, 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), see [19]. 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 a 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 high resolution 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 [1] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, New York, 1972, p. 11. [2] A. Collet, Separation and purification of enantiomers by crystallization methods, Enantiomer 4 (1999) 157–172. [3] M.P. Elsner, D. Fernández Menéndez, E. Alonso Muslera, A. Seidel-Morgenstern, Experimental study and simplified mathematical description of preferential crystallization, Chirality 17 (2005) 183–195. [4] M.P. Elsner, G. Ziomek, A. Seidel-Morgenstern, Simultaneous preferential crystallization in a coupled batch operation mode. Part I: Theoretical analysis and optimization, Chem. Eng. Sci. 62 (2007) 4760–4769. [5] R. Gunawan, I. Fusman, R.D. Braatz, High resolution algorithms for multidimensional population balance equations, AIChE J. 50 (2004) 2738–2749. [6] M.J. Hounslow, A discretized population balance for continuous systems at steady-state, AIChE J. 36 (1990) 106–116. [7] J. Jacques, A. Collet, S.H. Wilen, Enantiomers, Racemates and Resolutions, Wiley, New York, 1994. [8] B. Koren, A robust upwind discretization method for advection, diffusion and source terms, in: C.B. Vreugdenhill, B. Koren (Eds.), Numerical Methods for Advection–Diffusion Problems, Notes on Numerical Fluid Mechanics, vol. 45, Vieweg Verlag, Braunschweig, 1993, pp. 117–138 (Chapter 5). [9] S. Kumar, D. Ramkrishna, On the solution of population balance equations by discretization–III. Nucleation, growth and aggregation of particles, Chem. Eng. Sci. 52 (1997) 4659–4679. [10] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge Univ. Press, New York, NY, 2002. [11] Y.I. Lim, J.-M.L. Lann, X.M. Meyer, X. Joulia, G. Lee, E.S. Yoon, On the solution of population balance equation (PBE) with accurate front tracking method in practical crystallization processes, Chem. Eng. Sci. 57 (2002) 3715–3732. [12] A. Ma, D.K. Tafti, R.D. Braatz, High resolution simulation of multidimensional crystal growth, Ind. Eng. Chem. Res. 41 (2002) 6217–6223. [13] S. Qamar, A. Ashfaq, M.P. Elsner, I. Angelov, G. Warnecke, A. Seidel-Morgenstern, Adaptive high resolution schemes for multidimensional population balances in crystallization processes, Comp. & Chem. Eng 31 (2007) 1296–1311. [14] S. Qamar, M.P. Elsner, I. Angelov, G. Warnecke, A. Seidel-Morgenstern, A comparative study of high resolution schemes for solving population balances in crystallization, Comp. & Chem. Eng. 30 (2006) 1119–1131. [15] D. Ramkrishna, Statistical models of cell populations, Adv. Biochem. Engrg. 11 (1979) 1–47. [16] D. Ramkrishna, Population Balances: Theory and Applications to Particulate Systems in Engineering, Academic Press, San Diego, CA, 2000.
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
753
[17] R.A. Sheldon, L.S. Hulshof, A. Bruggink, F.J.J. Leusen, A.D. van der Haest, H. Wijnberg, Crystallization techniques for the industrial synthesis of pure enantiomers, in: Proceedings of the Chiral 90 Symposium, Manchester, United Kingdom, Spring Innovations Ltd., Stockport, England, 1990, pp. 101–107. [18] B. van Leer, Upwind-difference methods for aerodynamic problems governed by the Euler equations, in: Lectures in Applied Mathematics, vol. 22, American Mathematical Society, Providence, RI, 1985, pp. 327–336. [19] X.Y. Woo, R.B.H. Tan, P.S. Chow, R.D. Braatz, Simulation of mixing effects in antisolvent crystallization using a coupled CFD-PDF-PBE approach, Crystal Growth & Design 6 (2006) 1291–1303.
Numerical approximations of a population balance model for coupled batch preferential crystallizers S. Qamar a,∗ , I. Angelov b , M.P. Elsner b , A. Ashfaq a , A. Seidel-Morgenstern b , G. Warnecke a a Institute for Analysis and Numerics, Otto-von-Guericke University, PSF 4120, 39106 Magdeburg, Germany b Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstrasse 1, 39106 Magdeburg, Germany
Available online 1 April 2008
Abstract This article is concerned with the numerical approximations of population balance equations for modeling coupled batch preferential crystallization processes. The current setup consists of two batch crystallizers interconnected with two fines dissolution pipes. The crystallization of both enantiomers is assumed to take place in separate crystallizers after seeding with their corresponding crystals. The withdrawn fines are assumed to be dissolved in the dissolution unit after heating. crystallizer, the crystallizer temperature before entering to the opposite crystallizer. Two types of numerical methods are proposed for the simulation of this process. The first method uses high resolution finite volume schemes, while the second method is the so-called method of characteristics. 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 method of characteristics 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. A numerical test problem with both isothermal and non-isothermal conditions is considered here. The numerical results show clear advantages of the proposed schemes. © 2008 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Population balance models; Preferential crystallization of enantiomers; Coupled crystallizers; Fines dissolution; High resolution schemes; Method of characteristics; Nucleation and growth
1. Introduction Separation of chiral molecules is an important industrial process as many, especially organic, molecules are chiral. Usually, one of the enantiomers shows desired properties for therapeutic activities or metabolism, whereas the other enantiomer may be inactive or may even cause some undesired effects. An attractive and energy efficient way for separation of chiral substances is the enantioselective preferential crystallization. This concept is usually applied to the so-called conglomerates, i.e. a physical mixture of crystals where each crystal is enantiomerically pure. Such systems tend to reach an equilibrium state in which the liquid phase will have racemic composition and the solid * Corresponding author. Tel: +49 391 6712027; fax: +49 391 6718073.
E-mail address: [email protected] (S. Qamar). 0168-9274/$30.00 © 2008 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2008.03.033
740
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Fig. 1. Principle of preferential crystallization illustrated in ternary phase diagram.
phase will consist of a mixture of crystals of both enantiomers. However, before approaching this state, it is possible to preferentially produce just one of the enantiomers after seeding with homochiral crystals. The principle of a preferential crystallization process can be illustrated in a ternary-phase diagram. Fig. 1 depicts the solvent and two pure enantiomers in the corners of the diagram. Starting from a saturated solution at temperature Tcryst +T , the solution becomes supersaturated if it is rapidly cooled down to temperature Tcryst within the metastable zone. In this zone the liquid will remain free of particles for a limited amount of time, as no spontaneous (primary) nucleation will occur. In Fig. 1 point A represents the initial mixture of two enantiomers and solvent (e.g.: L-threonine, D-threonine and water). An enantiomeric excess (as depicted) is beneficial for the separation process, but not strictly necessary. If at point A the crystallizer is seeded with seeds of both enantiomers E1 and E2 , the two crystal types will start to grow, and induce secondary nucleation. The process will eventually end at point E, which represents the equilibrium point of the two enantiomers for the temperature Tcryst . However, this is not an example of enantioselective preferential crystallization, as the two enantiomers are being produced simultaneously. On the other hand, if at point A the crystallizer is seeded with homochiral seeds of type E1 only, the liquid phase composition tends towards point M. However, the experiments show that after certain time, spontaneous crystallization of the unseeded enantiomer (in this case E2 ) is observed, and therefore the trajectory is attracted again by the common equilibrium point E, see [3,4]. Therefore the crystallization conditions have to be designed carefully and the process must be stopped before a significant nucleation of the unseeded (counter) enantiomer occurs – at point B for example. Detailed treatises of the process of crystallization by entrainment can be found in the literature [2,7,17]. In this article, the crystallization of the two enantiomers is assumed to take place in two separate crystallizers, coupled by their fines dissolution loops, see Fig. 2. The extracted solution is screened by filters and assumed to be free of larger crystals. Therefore, only small particles are withdrawn to the fines dissolution loop. In order to assure a crystal-free liquid exchange, the withdrawn liquid in the fines dissolution loop is heated, so that the liquid becomes undersaturated and the withdrawn small particles dissolve. Before entering into the opposite crystallizer, this liquid is assumed to be cooled again. This ensures that there will be no negative effect on the particles inside the opposite crystallizer due to warm liquid flux. The attrition or agglomeration processes are not considered in the current study. Furthermore, it is assumed that crystals of the one enantiomer do not change the birth and nucleation rates of the counter enantiomer. Population balance models (PBMs) are the widely used simulation tool for different processes encountered in several scientific and engineering disciplines, see [15]. They can be used to describe the time evolution of one or more property distributions of individuals population. PBMs are used to study e.g., precipitation, polymerization, crystallization, particle size distribution (PSD) of crushed material and rain drops, dispersed phase size 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.
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
741
Fig. 2. Coupled batch process setup.
Since population balance equations (PBEs) can be solved analytically only for simplified cases, numerical schemes are usually needed. Many studies are therefore focused on the development of accurate and efficient numerical solutions of the population balance equations, see Ramkrishna [16]. High resolution schemes, which were originally developed for gas dynamics, were also applied for the numerical solution of PBEs, see Ma et al. [12], Gunawan et al. [5]. Qamar et al. [13,14] have also used high resolution finite volume schemes in order to solve one and two-dimensional population balance equations 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 [9] 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. [11] applied high resolution spatial discretization methods (WENO schemes) and the method of characteristics 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 coupled batch preferential crystallization process. The model is further elaborated by considering the isothermal and non-isothermal conditions. We have implemented the high resolution scheme of Koren [8] and that by LeVeque [10] for solving PBEs which simulate preferential crystallization process in two coupled crystallizers. Furthermore, we have also implemented the method of characteristics [9,11]. 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 RK4 method is used, which is a Runge–Kutta method of order four. The paper is organized as follows. In Section 2, we start with the model for two coupled batch preferential crystallizers which are interconnected with fines dissolution units. In Section 3, we briefly derive our proposed numerical schemes for solving the current PBEs. Section 4 includes a numerical test problem for the current model for isothermal and non-isothermal cases. In Section 5, we give conclusions and remarks. 2. Coupled preferential crystallizers model with fines dissolution In this section, we provide a mathematical model for the simulation of two batch preferential crystallizers which are interconnected with fines dissolution pipes. For more details about preferential crystallization the reader is referred to [3,4]. A simplified dynamic model of ideally mixed batch crystallizers for isothermal and non-isothermal conditions is simulated. Here, we are considering the simplest case, in which inflow and outflow pipes are attached to the crystallizers and we assume that fines dissolves completely at the end of the pipe. Furthermore, we are assuming the same volume and residence time for both crystallizers. Similar assumptions hold for the two dissolution pipes as their
742
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Fig. 3. Death function h(x).
volume and residence time are identical. The function of pipes (heat exchanger) is to provide sufficient heat to the fines so that they can dissolve completely. Fig. 2 shows the schematic diagram of the coupled preferential crystallization process. In vessel A the L-threonine is the preferred enantiomer, D-threonine is the counter enantiomer and in vessel B vice versa. The population balance equations describing the balance laws for the solid phase are as follows: (k)
(k)
1 ∂nα (x, t) ∂nα (x, t) = −G(k) − hα (x)n(k) α (t) α (x, t) ∂t ∂x τ1,α with initial data
(2.1)
(k)
n(k) α (x, 0) = n0,α (x),
(2.2)
where k ∈ {p, c} and α ∈ {A, B}. The notation k stands for preferred or counter enantiomer and α for crystallizer A or B. In Eq. (2.1) τ1,α is the residence time in the corresponding crystallizer. It is defined as Vα , (2.3) τ1,α = V˙α where Vα is the volume of each crystallizer and V˙α is the corresponding volumetric flux rate. In Eq. (2.1) hα (x) is a death function describing the dissolution of particles below some critical size. It can be defined as a step (Heaviside) function. Fig. 3 represents one possible death function which we have considered in our numerical test problems with xcrit = 2 × 10−4 m. The mass balances for the liquid phase in the crystallizers can be described by the following equations (k)
dmA (t) (k) =m ˙ (k) ˙ (k) out,A (t) − 3ρkv GA (t) in,A (t) − m dt
∞
x 2 n(k) A (x, t) dx,
(2.4)
0 (k) dmB (t)
dt
(k)
(k)
(k)
∞
=m ˙ in,B (t) − m ˙ out,B (t) − 3ρkv GB (t)
(k)
x 2 nB (x, t) dx
(2.5)
0
with initial data (k)
(k)
mA (0) = m0,A ,
(k)
(k)
mB (0) = m0,B .
(2.6)
Here, ρ is the crystal density and kv is the volume shape factor defined such that the volume of a crystal with length (k) x is kv x 3 . Because of the fines dissolution Eqs. (2.4) and (2.5) have four mass fluxes. The first one m ˙ in,A (t) is the (k)
incoming mass flux in the liquid phase into the crystallizer A. The second one m ˙ out,A (t) is the outgoing mass flux from (k)
(k)
crystallizer A into the dissolution pipe. The flux m ˙ in,B (t) is the incoming mass flux into crystallizer B and m ˙ out,B (t) is the outgoing flux from crystallizer B. They are defined as follows (k) ˙ m ˙ (k) out,α (t) = wα ρliq Vα ,
k ∈ {p, c}, α ∈ {A, B}, ∞ kv ρ (k) (k) m ˙ in,α (t) = m ˙ out,α (t − τ2 ) + x 3 hα (x)n(k) α (x, t − τ2 ) dx, τ1,α 0
(2.7) (2.8)
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
743
(k)
where wα are the mass fractions and ρliq is the liquid density. Here τ2 is the residence time of the dissolution unit (the pipe) interconnecting both crystallizers. It is defined as τ2 =
πr 2 l , V˙
(2.9)
where r and l are the radius and length of each pipe, respectively. Due to the time delay in (2.8), special attention is needed for the numerical approximation of this equation. To simplify our problem, τ2 is chosen to be an integer multiple of the fixed time step t. At each time step the values of (k) the mass flow rates to the dissolution pipes mout,α (t) are saved in vectors which can be then used in (2.8) for t − τ2 0. (k) Similarly, The values of nα (x, t) are saved in matrices for using in (2.8) for t − τ2 0. For the case t − τ2 < 0 we (k) take min,α (t) = 0. However, one can also use an adaptive ODE-solver and some interpolation formula to calculate the terms in (2.8). It should be stressed that at the beginning of the process, i.e. t = 0, the two tanks including pipes are treated as one unit. For this reason the initial mass fractions for both enantiomers are everywhere the same. The birth of new particles (nucleation) is incorporated in this model as a left boundary condition for the PDE at minimum crystal size x = x0 (k)
n(k) α (x0 , t) =
Bα (t) (k)
(2.10)
.
Gα (t)
The assumption for this substitution is that if Bα is the rate of appearance of nearly zero-sized particles and the growth rate is independent of size, then dNα dNα dx (k) · = = nα (x0 , t)G(k) (2.11) Bα (t) = α (t). dt x→x0 dx dt x→x0 Here Nα is the number of crystals formed in crystallizer α ∈ {A, B}. The growth rate kinetics are described as: (k) g G(k) α (t) = kg Sα (t) − 1 .
(2.12) (k)
Here kg is the growth rate constant and g is the growth rate exponent. The function Sα denotes the relative supersaturation for the corresponding enantiomer defined as: (k)
Sα(k) (t) =
wα (t) (k)
(2.13)
.
weq,α (t)
(k)
(k)
The term wα (t) is mass fraction of the dissolved enantiomer, while weq,α (t) is equilibrium mass fraction for the corresponding enantiomer. The mass fractions are defined as: (k)
wα(k) =
mα (p)
(c)
mα + mα + mW
(2.14)
,
where mW is the mass of solvent (here water). The nucleation rate kinetics are described as: b(p) (p) (p) (p) μ3,α (t), Sα (t) − 1
Bα(p) (t) = kb
Bα(c) (t) = kb(c) e
−
b(c) (c) ln(Sα (t))2
,
(2.15)
where kb(k) are the nucleation rate constants and b(k) are the nucleation rate exponents for k ∈ {p, c}. Here, μ(k) 3,α (t) (k)
represent the third moments of the number densities nα (x, t). The given model covers the case of isothermal operation, i.e. the temperature is kept constant during the batch process. A natural extension of the process would be to vary the temperature during the batch process according to a specified temperature profile.
744
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Fig. 4. Cell centered finite volume grid.
3. Numerical schemes In this section we introduce two types of numerical methods for the solution of the above population balance models 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 ∂[G(x, t)n(x, t)] ∂n(x, t) =− + B0 (t) δ(x − x0 ) + Q(x, t), (3.1) ∂t ∂x 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 [6]. 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 homogeneous PBE and defining the ratio of nucleation and growth terms as a left boundary condition, i.e. ∂[G(x, t)n(x, t)] ∂n(x, t) =− + Q(x, t), ∂t ∂x
with n(x0 , t) =
B0 (t) . G(t)
(3.2)
3.1. High resolution schemes Here, we give the formulation of high resolution finite volume schemes for the numerical solutions of the above PBE. For a discussion of the schemes the reader is referred to our article [14]. 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 denote by (xi− 1 )i∈{1,...,N +1} a mesh of [x0 , xmax ], where x0 2 is the minimum and xmax is the maximum crystal length of interest. As shown in Fig. 4, for each i = 1, 2, . . . , N , x represents the cell width, the points xi refer to the cell centers, and the points xi± 1 represent the cell boundaries. 2 The integration of (3.2) over the cell Ωi = [xi− 1 , xi+ 1 ] yields the following cell centered semidiscrete finite2 2 volume scheme for Fi± 1 := (Gn)i± 1 2 2 ∂n(x, t) dx = −(Fi+ 1 − Fi− 1 ) + Q(x, t) dx. (3.3) 2 2 ∂t Ωi
Ωi
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.4) ni (t) = x x Ωi
Ωi
Using Eq. (3.4), Eq. (3.3) implies 1 dni =− (F 1 − Fi− 1 ) + Qi , for i = 1, 2, . . . , N (3.5) 2 dt x i+ 2 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.
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Fx
i+ 21
= Fi ,
Fx
i− 21
= Fi−1 .
745
(3.6)
High order accuracy is easily obtained by piecewise polynomial interpolation. One can take for instance 1+κ 1−κ (3.7) (Fi+1 − Fi ) + (Fi − Fi−1 ), 4 4 where κ ∈ [−1, 1]. Similarly one can write an expression for Fi− 1 . Here κ is a parameter that has to be chosen from 2 the 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 by Koren [8] and uses the value κ = 13 , while the second scheme was given by LeVeque [10] 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 [8,14]. Fi+ 1 = Fi + 2
HR-κ = 1/3 scheme: According to this scheme the flux at the right boundary xi+ 1 of the cell Ωi in (3.5) is [8] 2
1 Fi+ 1 = Fi + Φ(ri+ 1 )(Fi − Fi−1 ), 2 2 2 where the flux limiting function Φ according to Koren [8] is defined as 1 2 + r 1,2 . Φ(ri+ 1 ) = max 0, min 2ri+ 1 , min 2 2 3 3 i+ 2
(3.8)
(3.9)
Here the argument ri+ 1 of this function is the so-called upwind ratio of two consecutive flux gradients 2
ri+ 1 = 2
Fi+1 − Fi + ε . Fi − Fi−1 + ε
(3.10)
This expression has to be evaluated with a small parameter, e.g. ε = 10−10 , to avoid division by zero [8]. 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) =
B0 (t) . Gx −
(3.11)
1
HR-κ = −1 scheme: In this scheme the expression for the numerical flux can be obtained from Eq. (3.7) by taking κ = −1. The flux at the right boundary xi+ 1 of the cell Ωi in (3.5) is [10] 2
1 Fi+ 1 = Fi + φ(θi+ 1 )(Fi+1 − Fi ), 2 2 2
where θi+ 1 = 2
Fi − Fi−1 . Fi+1 − Fi
(3.12)
Here the limiting function φ is given by the van Leer flux limiter [18] φ(θi+ 1 ) = 2
|θi+ 1 | + θi+ 1 2
2
1 + |θi+ 1 |
(3.13)
2
with the same boundary condition as given by (3.11). 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 scheme (3.6) 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.6)–(3.12) have to be mirrored at the considered volume boundaries xi+ 1 for i = 1, 2, . . . , N . Finally, the resultant 2 system of ordinary differential equations can be solved by any standard ODEs solver. In the current study we have used the RK4 method which is a Runge–Kutta method of order four. For approximating the integral terms in the mass balance equations (2.4) and (2.5) we have used the composite Simpson’s rule.
746
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
3.2. Method of characteristics (MOC) For linear advection equations a good alternative to the above schemes is the method of characteristics. 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 points xi+1/2 (t) for i = 0, 1, 2, . . . , N . Let us substitute the growth rate G(t, x) by dx := G(x, t). dt Then Eq. (3.2) leads to ∂n(x, t) ∂ dx + · n(x, t) = Q(x, t). ∂t ∂x dt
(3.14)
(3.15)
Integration over the control volume Ωi (t) = [xi− 1 (t), xi+ 1 (t)] gives 2
Ωi (t)
2
x 1 (t) i+ dx ∂n(x, t) dx + · n(x, t) 2 = Q(x, t) dx. ∂t dt x 1 (t) i− 2
Ωi (t)
The Leibniz formula [1] on the left-hand side of (3.16) gives d n(x, t) dx = Q(x, t) dx. dt Ωi (t)
(3.16)
(3.17)
Ωi (t)
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.4), they are defined as 1 1 n(x, t) dx, Qi := Q(x, t) dx, (3.18) ni := xi (t) xi (t) Ωi
Ωi
where xi (t) = xi+ 1 (t) − xi− 1 (t). After using (3.18), Eq. (3.17) gives 2
2
d xi+ 1 (t) − xi− 1 (t) ni = xi (t) Qi . (3.19) 2 2 dt By using the product rule and (3.14), the left-hand side of (3.19) gives dx 1 dxi− 1 dni d i+ 2 2 xi+ 1 (t) − xi− 1 (t) ni = xi (t) + − ni 2 2 dt dt dt dt dni + (Gi+ 1 − Gi− 1 )ni . = xi (t) (3.20) 2 2 dt After replacing the left-hand side of (3.19) with (3.20) and dividing the resultant equation by xi (t) one gets dni ni = −(Gi+ 1 − Gi− 1 ) + Qi . 2 2 dt xi (t) In summary to use MOC, we have to solve the following set of equations dni ni = −(Gi+ 1 − Gi− 1 ) + Qi , 2 2 xi (t) dt dxi+ 1 2 = Gi+ 1 , ∀i = 1, 2, . . . , N 2 dt
(3.21)
(3.22) (3.23)
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
747
with initial data n(xi , 0) = n0 (xi ). Note that, in case of size independent growth the first term on the right-hand side of (3.22) 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 the case of stiff nucleation a small time step or comparatively refined mesh in the vicinity of the nucleation is desirable to avoid the loss in mass and unphysical oscillations, see [11]. For using different mesh speeds near the left boundary one can use the technique introduced in [9]. Usually, there is no effect on the MOC when the mesh becomes more refined, except an increase of the computational time. Here, the boundary condition is again the same as given by (3.11). The above system of ordinary differential equations can be solved by any standard ODE solver. In the current study we have used the RK4 method which is a Runge–Kutta method of order four. For approximating the integral terms in the mass balance equations (2.4) and (2.5) we have used the composite Simpson’s rule. 4. Numerical test problems In order to validate our numerical schemes for the current model for preferential crystallization, we consider the following numerical test problem. In this problem we have taken the minimum crystal size x0 = 0. The initial number density of the preferred enantiomer in each crystallizer is given as 1 ln(x) − μ 2 kV ρ (p) (p) , Iα = · exp − μ (0). (4.1) nα (x, 0) = √ √ Ms,α 3,α x 2π σ Iα 2σ The corresponding number density of the counter enantiomer in each crystallizer is zero, i.e. n(c) α (x, 0) = 0.
(4.2)
We assume σ = 0.3947 m, μ = −6.8263, and Ms,α are the masses of the initial seeds. The maximum crystal size that was expected is xmax = 0.005 m which is subdivided into 500 grid points. The final time for the simulation was taken as 600 minutes. The kinetic parameters considered in this problem are given in Table 1. They are taken to describe the crystallization of the enantiomers of amino acid threonine in water. Note that enantiomer E1 is preferred in crystallizer A and enantiomer E2 is preferred in crystallizer B. The temperature trajectory used for the simulation is as follows T (t)[C o ] = −1.24074 × 10−7 t 3 + 4.50926 × 10−5 t 2 − 0.00405556t + 33. Here the time is taken in minutes. For this temperature trajectory, the constants kg and temperature as given below EA,g
kg = kg,0 e− R(T +273.15) ,
EA,b
kb = kb,0 e− R(T +273.15) . (k)
(k)
(4.3) (p) kb
will become functions of
(4.4)
(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 the temperature. They are given as (p)
weq,α (t) =
2
T i (t) Ai + Bi w (c) (t) ,
(4.5)
T i (t) Ci + Di w (p) (t) .
(4.6)
i=0 (c) weq,α (t) =
2 i=0
The constants in the above two equations are given in Table 1. Note that the same Eqs. (4.5) and (4.6) are also used in the isothermal case (the temperature remains constant). Fig. 5 shows the number density plots in tank A for the preferred enantiomer with isothermal 33 ◦ C and for the non-isothermal condition, temperature follows Eq. (4.3). Note that, E1 is the preferred enantiomer in Tank A and E2 is the preferred enantiomer in Tank B. Due to the same initial conditions (mass of seeds, initial CSD) and due to the phenomenon that enantiomers have identical properties in arhiral environments, the plots for the enantiomer E1 in
748
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Table 1 Parameters for test problem 1 and 2 Description
Symbols
Value
Unit
Growth rate constant Growth rate exponent Activation energy
kg,0 g EA,g
4.62 × 108 1.0 75.6 × 103
m/s – kJ/mol
Nucleation rate constant (seeded)
kb,0
(p)
3.24 × 1025
1/m3 s
Nucleation rate exponent (seeded) Activation energy
b(p)
4.0 78.7 × 103
– kJ/mol
Nucleation rate constant (unseeded)
kb,0
(c)
3.84 × 106
1/s
Nucleation rate exponent (unseeded) Universal gas constant Density of crystals Density of water Volume shape factor Initial mass of p-enantiomer
b(c)
6.0 × 10−2 8.314 1250 1000 0.0248 0.100224
– J/mol K kg/m3 kg/m3 – kg
Initial mass of c-enantiomer Mass of water Mass of seeds Residence time in crystallizer Volume of crystallizer Residence time 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)
m0,α mW Ms τ1 V1 τ2 r l A0 A1 A2 B0 B1 B2 C0 C1 C2 D0 D1 D2
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
kg kg kg min m3 min m m – – – – – – – – – – – –
EA,b
R ρ ρw kv (p) m0,α (c)
Table 2 Mass preservation in the schemes Method
Percentage error isothermal case
Percentage error non-isothermal case
First order scheme HR-κ = −1 scheme HR-κ = 1/3 scheme MOC
2.945 3.020 3.022 1.781
2.148 2.301 2.302 1.996
tank A and for the enantiomer E2 in tank B are the same. In the isothermal case, the first large peak in the number density distribution is due to high nucleation rate of the preferred enantiomer at the beginning of the process. While the second small number density peak is due to the initial crystals growth. In the non-isothermal case, the temperature is a decreasing function of time. Hence the peak in the number density generated by nucleation has a higher value as compared to the isothermal case because of the temperature profile. The number of peaks in the crystal size distribution for the non-isothermal case entirely depends on the temperature profile used. The small zigzags in the left-hand side plot of MOC is due to stiff nucleation which was treated as a boundary condition. In MOC the mesh moves to the right at each time step and we have to add one mesh cell at the left boundary to include the boundary conditions. In the case
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
749
Fig. 5. Comparison of the number density for the preferred enantiomer for different operation modes at t = 600 minutes.
Fig. 6. Three dimensional plots of number density for the preferred enantiomer by using MOC.
of stiff nucleation a small time step or comparatively refined mesh in the vicinity of the nucleation is desired to avoid such oscillations, see [9,11]. There is no visible difference between results of the given schemes, so one cannot say exactly which scheme is better one. Table 2 shows that percentage errors for mass preservation are minimum for the MOC as compared to the finite volume schemes. The reason can be the current stiff nucleation which is considered as a boundary condition in the finite volume schemes. In the case of high resolution schemes this stiff nucleation may produce over predictions. This may lead to the bigger errors in the mass balances. Furthermore, there is also no significant difference in mass fractions of the MOC and the finite volume schemes, as well as supersaturation, growth rate, nucleation rates and the third moment. Therefore we have only presented here the plots of the MOC. Figs. 6 and 7 show the three dimensional plot of number density for the preferred enantiomer and the counter enantiomer. Fig. 8 shows the mass fraction plots for the preferred (p-) and counter (c-) enantiomers. In the isothermal case, the mass fraction of the p-enantiomer decreases sharply because we have seeded the p-enantiomer and it crystallizes out. While c-enantiomer mass fraction stay 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. 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.5) and (4.6). Figs. 10 and 11 show the growth rate and nucleation rate plots. At the beginning, growth rate and nucleation rate for the p-enantiomer are reduced significantly because of sharp changes in mass fraction and supersaturation. In the non-isothermal case, nucleation rate for counter enantiomer achieves a very high value because of the high supersaturation level.
750
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Fig. 7. Three dimensional plots of number density for the counter enantiomer using MOC.
Fig. 8. Comparison of the mass fractions for the preferred (p-) and counter (c-) enantiomer using MOC.
Fig. 9. Comparison of the supersaturation for the preferred (p-) and counter (c-) enantiomer using MOC.
Fig. 12 shows the trajectories of the third moments. These plots in the isothermal case for the p-enantiomer stays constant at the end of the process, because of no further change in mass fraction and supersaturation. In the nonisothermal case the third moment does not approach a steady state because of the decreasing temperature in the crystallizer.
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
Fig. 10. Comparison of growth rate for the preferred (p-) and counter (c-) enantiomer using MOC.
Fig. 11. Comparison of nucleation rate for the preferred (p-) and counter (c-) enantiomer using MOC.
Fig. 12. Comparison of the third moment for the preferred (p-) and counter (c-) enantiomer using MOC.
751
752
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
5. Conclusions In this article two high resolution finite volume schemes and a method of characteristics were proposed for solving population balance models in two coupled preferential crystallizers. For the process of preferential crystallization there are two main advantages of considering two coupled crystallizers which are interconnected by two fines dissolution units. The first one is that one gets both enantiomers at the same time in separate crystallizers. Secondly, because of the fines dissolution, the amount of small particles reduces which further enhances the particle growth. For handling the nucleation term in the method of characteristics, a cell of nuclei size is added at a given time level. The resultant system of ODEs are then solved by a standard ODE-solver. In this work we have used the RK4 method which is a Runge–Kutta method of order four. Both high resolution schemes and the method of characteristics worked very well for the numerical test problems which are considered in this article. However, it was found that the method of characteristics resolves better the number density distribution and give smaller errors 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 method of characteristics. The use of MOC for the current growth processes looks quite successful. However, due to their generality and robustness to model perturbations, 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), see [19]. 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 a 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 high resolution 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 [1] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, New York, 1972, p. 11. [2] A. Collet, Separation and purification of enantiomers by crystallization methods, Enantiomer 4 (1999) 157–172. [3] M.P. Elsner, D. Fernández Menéndez, E. Alonso Muslera, A. Seidel-Morgenstern, Experimental study and simplified mathematical description of preferential crystallization, Chirality 17 (2005) 183–195. [4] M.P. Elsner, G. Ziomek, A. Seidel-Morgenstern, Simultaneous preferential crystallization in a coupled batch operation mode. Part I: Theoretical analysis and optimization, Chem. Eng. Sci. 62 (2007) 4760–4769. [5] R. Gunawan, I. Fusman, R.D. Braatz, High resolution algorithms for multidimensional population balance equations, AIChE J. 50 (2004) 2738–2749. [6] M.J. Hounslow, A discretized population balance for continuous systems at steady-state, AIChE J. 36 (1990) 106–116. [7] J. Jacques, A. Collet, S.H. Wilen, Enantiomers, Racemates and Resolutions, Wiley, New York, 1994. [8] B. Koren, A robust upwind discretization method for advection, diffusion and source terms, in: C.B. Vreugdenhill, B. Koren (Eds.), Numerical Methods for Advection–Diffusion Problems, Notes on Numerical Fluid Mechanics, vol. 45, Vieweg Verlag, Braunschweig, 1993, pp. 117–138 (Chapter 5). [9] S. Kumar, D. Ramkrishna, On the solution of population balance equations by discretization–III. Nucleation, growth and aggregation of particles, Chem. Eng. Sci. 52 (1997) 4659–4679. [10] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge Univ. Press, New York, NY, 2002. [11] Y.I. Lim, J.-M.L. Lann, X.M. Meyer, X. Joulia, G. Lee, E.S. Yoon, On the solution of population balance equation (PBE) with accurate front tracking method in practical crystallization processes, Chem. Eng. Sci. 57 (2002) 3715–3732. [12] A. Ma, D.K. Tafti, R.D. Braatz, High resolution simulation of multidimensional crystal growth, Ind. Eng. Chem. Res. 41 (2002) 6217–6223. [13] S. Qamar, A. Ashfaq, M.P. Elsner, I. Angelov, G. Warnecke, A. Seidel-Morgenstern, Adaptive high resolution schemes for multidimensional population balances in crystallization processes, Comp. & Chem. Eng 31 (2007) 1296–1311. [14] S. Qamar, M.P. Elsner, I. Angelov, G. Warnecke, A. Seidel-Morgenstern, A comparative study of high resolution schemes for solving population balances in crystallization, Comp. & Chem. Eng. 30 (2006) 1119–1131. [15] D. Ramkrishna, Statistical models of cell populations, Adv. Biochem. Engrg. 11 (1979) 1–47. [16] D. Ramkrishna, Population Balances: Theory and Applications to Particulate Systems in Engineering, Academic Press, San Diego, CA, 2000.
S. Qamar et al. / Applied Numerical Mathematics 59 (2009) 739–753
753
[17] R.A. Sheldon, L.S. Hulshof, A. Bruggink, F.J.J. Leusen, A.D. van der Haest, H. Wijnberg, Crystallization techniques for the industrial synthesis of pure enantiomers, in: Proceedings of the Chiral 90 Symposium, Manchester, United Kingdom, Spring Innovations Ltd., Stockport, England, 1990, pp. 101–107. [18] B. van Leer, Upwind-difference methods for aerodynamic problems governed by the Euler equations, in: Lectures in Applied Mathematics, vol. 22, American Mathematical Society, Providence, RI, 1985, pp. 327–336. [19] X.Y. Woo, R.B.H. Tan, P.S. Chow, R.D. Braatz, Simulation of mixing effects in antisolvent crystallization using a coupled CFD-PDF-PBE approach, Crystal Growth & Design 6 (2006) 1291–1303.