Numerical solution of population balance ... - Semantic Scholar
Computers and Chemical Engineering 31 (2007) 1576–1589
Numerical solution of population balance equations for nucleation, growth and aggregation processes Shamsul Qamar ∗ , Gerald Warnecke Institute for Analysis and Numerics, Otto-von-Guericke University, PSF 4120, 39106 Magdeburg, Germany Received 26 August 2006; received in revised form 16 January 2007; accepted 16 January 2007 Available online 28 January 2007
Abstract This article focuses on the derivation of numerical schemes for solving population balance models (PBMs) with simultaneous nucleation, growth and aggregation processes. Two numerical methods are proposed for this purpose. The first method combines a method of characteristics (MOC) for growth process with a finite volume scheme (FVS) for aggregation process. For handling nucleation terms, a cell of nuclei size is added at a given time level. The second method purely uses a semi-discrete finite volume scheme for nucleation, growth and aggregation of particles. Note that both schemes use the same finite volume scheme for aggregation process. On one hand, the method of characteristics offers a technique which is in general a powerful tool for solving linear growth processes, has the capability to overcome numerical diffusion and dispersion, is computationally efficient, as well as give highly resolved solutions. On the other hand, the finite volume schemes which were derived for a general system in divergence form, are applicable to any grid to control resolution, and are also computationally not expensive. In the first method a combination of finite volume scheme and the method of characteristics gives a highly accurate and efficient scheme for simultaneous nucleation, growth and aggregation processes. The second method demonstrates the applicability, generality, robustness and efficiency of high-resolution schemes. The proposed techniques are tested for pure growth, simultaneous growth and aggregation, nucleation and growth, as well as simultaneous nucleation, growth and aggregation processes. The numerical results of both schemes are compared with each other and are also validated against available analytical solutions. The numerical results of the schemes are in good agreement with the analytical solutions. © 2007 Elsevier Ltd. All rights reserved. Keywords: Population balance models; High-resolution finite volume schemes; Method of characteristics; Nucleation; Growth; Aggregation; Hyperbolic conservation laws; Particles
1. Introduction The dynamic behavior of a population of small particles is a subject of interest in various engineering fields. These areas include atmospheric physics, precipitation, crystallization, pharmaceutical manufacture, aerosol formation, colloid chemistry, growth of microbial and cell populations, and so on. Most of such systems involve simultaneous nucleation, growth and aggregation of particles. Population balances, in their integrodifferential equation form, are a widely used tool for simulating these processes. Since, population balance equations (PBEs) can only be solved exactly in simplified cases, numerical solutions are usually needed. Many studies are therefore focused on
∗
Corresponding author. Tel.: +49 391 6712027; fax: +49 391 6718073. E-mail address: [email protected] (S. Qamar). 0098-1354/$ – see front matter © 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2007.01.006
the accurate and efficient numerical solutions of the population balance equations. In the past several numerical techniques were introduced for solving PBEs, see for example Ramkrishna (1985) for a review. Discussion on the stability and numerical dispersion of various finite difference approximations can be found in Lapidus and Pinder (1982). An extension of the finite difference-type discretization methods have been proposed by Hounslow, Ryall, and Marshall (1998) and David, Villermaux, and Klein (1991). However, their numerical results were found to be unsatisfactory, see Kumar and Ramkrishna (1997). A second-order spatial discretization for the growth term as well as an accurate and complex discretization for the aggregation term on a geometric grid was proposed by Litster, Smith, and Hounslow (1994). Later on, Muhr, David, Villermaux, and Jezequel (1996) used a firstorder upwind scheme for the spatial discretization of the growth term. However, due to numerical diffusion a loss in accuracy near steep fronts or discontinuities was observed.
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
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 term in order to solve PBEs with simultaneous nucleation, growth and aggregation processes. They found that, in contrast to other numerical techniques, MOC avoids the numerical dissipation error caused by the growth term discretization. For handling the nucleation term, which is usually difficult to treat, MOC was combined with a procedure of adding a cell of the nuclei size at a given time level. A standard ordinary differential equation (ODE) solver was used to solve the resultant ODEs. A short overview of the methods of Hounslow et al. (1988) and David et al. (1991) is also given in their article. The finite volume schemes, which were originally developed for the gas dynamics, were also applied for the numerical solution of population balance models in crystallization process, see Gunawan, Fusman, and Braatz (2004) and Ma, Tafti, and Braatz (2002a, 2002b) and Motz, Mitrovi´c, and Gilles (2002). Further, Lim, Lann, and Meyer (2002) applied highresolution spatial discretization methods (WENO schemes) and the modified method of characteristics (MOC) for the dynamic simulations of the batch crystallization including nucleation, growth, aggregation and breakage kinetics. They observed that, steep moving fronts or discontinuities appearing in the solution of the crystallization models were captured well by both methods without any spurious oscillations. Apart from these methods, Filbet and Laurenc¸ot (2004) proposed a numerical scheme for single-component aggregation which is based on a conservative finite volume formulation. The authors have rewritten the population balance equation for aggregation problems in a form which can be readily solved by a finite volume scheme. This special reformulation was a great achievement which enables one to apply the finite volume scheme in the aggregation case as well. Later on, Qamar and Warnecke (2006) have extended the scheme in order to simulate two-component aggregation problems. Apart from that, Qamar, Elsner, and Angelov (2006) and Qamar et al. (2007) have applied high-resolution finite volume schemes for solving population balance models in crystallization processes, where they considered growth and nucleation kinetics while neglecting the aggregation and breakage processes. The aim of the current work is to numerically solve PBEs for simultaneous nucleation, growth and aggregation processes. Two methods are proposed in this article. The first method combines a method of characteristics for growth term (Kumar & Ramkrishna, 1997) and a finite volume scheme (Filbet & Laurenc¸ot, 2004; Qamar & Ashfaq, 2006) for aggregation term. In the second method, a semi-discrete finite volume scheme (FVS) is used for all processes, see Qamar et al. (2006, 2007) and Qamar and Warnecke (2006) for a review. Note that, in both methods the aggregation kinetics is solved with the same finite volume formulation. The main difference in both methods is the way they calculate the growth term, for example the first method uses MOC for the growth term and the advection term simply disappears from the governing equation, while the second method uses a finite volume scheme to discretize the advection term. The efficiency and accuracy of the resultant methods are analyzed by comparing their numerical results with each other and with avail-
1577
able analytical solutions. For handling the nucleation term in the first scheme, MOC is combined with a procedure of adding a cell of the nuclei size at each time level. The same procedure was also used in Kumar and Ramkrishna (1997). A standard ODEs solver can be used to solve the resultant ODEs. There are two main differences between the first method and the one used in Kumar and Ramkrishna (1997). In our current method the aggregation term is treated with a finite volume scheme, while in Kumar and Ramkrishna (1997) the authors have used the fixed pivot technique for this purpose. Secondly, our numerical method uses a reformulated PBE instead of the original one. This reformulation is necessary in order to apply the finite volume scheme. In case of reformulated PBE, our numerical method calculates the volume (mass) density instead of the number density. However, one can easily recover back the number density at the end of simulation. Similarly Lim et al. (2002) used the original PBE and a WENO scheme which is far different from the current finite volume schemes. Note 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 method of characteristics to deal with the hyperbolic nature of the growth term, are needed when accurate population balance model solutions are needed to compute real-time 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 particle size distribution can be found in the recent articles of Shi, El-Farra, and Mhaskar (2006) and Shi, Mhaskar, and El-Farra (2005). Although we are not considering control problems in this work, the current schemes can be used in this direction as well. The paper is organized as follows. In Section 2, we start with an integro-differential PBE and reformulate it to a form which can be readily used for the derivation of our proposed numerical methods. In Section 3, we give two different grid discretizations and then derive both numerical methods for the solution of PBEs for simultaneous nucleation, growth and aggregation processes. In Section 4, we give some numerical results for different combinations of the nucleation, growth and aggregation processes. For validation of our numerical results we use the available analytical solutions. In Section 5, we give conclusions and remarks. 2. Mathematical model In this section we will reformulate the underlying onedimensional PBE in a form which can be readily used to derive the proposed numerical methods. Consider a PBE of the form: ∂f (t, x) ∂[G(t, x)f (t, x)] + ∂t ∂x = Qagg (t, x) + S(t, x),
(t, x) ∈ R2+ ,
(2.1)
1578
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
f (0, x) = f0 (x),
x ∈ R+ ,
(2.2)
where R+ :=]0, +∞[, f (t, x) is the number density; G(t, x) the growth term; S(t, x) the rate at which particle of size x are nucleated; Qagg (t, x) is the aggregation reaction term which is defined as x 1 Qagg (t, x) = β(t, x − x , x ) f (t, x − x ) f (t, x ) dx 2 0
∞ −
β(t, x, x ) f (t, x) f (t, x ) dx .
xi f (t, x) dx,
i = 0, 1, 2, . . . ,
(2.4)
0
where the zeroth moment M0 (t) and the first moment M1 (t) represent the total number of particles and the total volume of particles, respectively. We want to consider the mass balance. Multiplying both sides of Eqs. (2.1)–(2.3) by x, we get x
∂f (t, x) ∂[G(t, x)f (t, x)] +x ∂t ∂x = xQagg (t, x) + xS(t, x), x ∈ R+ ,
xf (0, x) = xf0 (x),
xQagg (t, x) =
1 2
x
(t, x) ∈ R2+ ,
(2.5) (2.6)
xβ(t, x − x , x ) f (t, x − x ) f (t, x ) dx
0
∞ −
xβ(t, x, x ) f (t, x) f (t, x ) dx .
(2.9)
where x ∞ u β(t, u, v) f (t, u) f (t, v) dv du,
F(t, x) = −
x ∈ R+ . (2.10)
˜ x) := xS(t, x). After Let us define f˜ (t, x) := xf (t, x) and S(t, using Eqs. (2.4)–(2.10) in (2.5) and (2.6), we obtain the mass balance: ∂f˜ (t, x) ∂[G(t, x)f˜ (t, x)] G(t, x)f˜ (t, x) + − ∂t ∂x x ∂F(t, x) ˜ =− + S(t, x), ∂x
f˜ (0, x) = f˜ 0 ,
(x, x ) ∈ R2+ .
It is basically the probability that particles of sizes x, x meet and stick together. The moments of the distribution function f are defined as Mi (t) =
∂F(t, x) , ∂x
0 x−u
The first term in Eq. (2.5) accounts for the formation of particles with volume x resulting from the merging of two particles with respective volume x and x − x , where x ∈ ]0, x[. The second term is the death term due to loss of particle of volumes x by aggregation with other particles of any size. The aggregation coefficient β = β(t, x, x ) characterizes the rate at which the aggregation of two particles with respective volumes x and x produces a particle of volume x + x and is a non-negative symmetric function,
∞
xQagg (t, x) = −
(2.3)
0
0 ≤ β(t, x, x ) = β(t, x , x),
Furthermore, the aggregation term xQagg (t, x) can be rewritten as (Filbet & Laurenc¸ot, 2004; Qamar & Ashfaq, 2006):
(2.7)
x ∈ R+ .
(2.12)
Instead of the original equation (2.1) we will discretize (2.11). Due to the relation f˜ (t, x) := xf (t, x) one can easily recover f (t, x). 3. Numerical procedure In the following we will derive two methods in order to solve the above modified population balance equations (2.11) and (2.12) for the simulation of simultaneous nucleation growth and aggregation processes. 3.1. Domain discretization In order to apply any numerical scheme, the first step is to discretize the computational domain. Here, we give two types of grid discretization, that is, regular (or irregular grid) and geometric grid. 3.1.1. Regular/irregular grid Let N be a large integer representing the total number of cells, we denote by (xi−(1/2) )i ∈ {1,...,N+1} a mesh of [xmin , xmax ], where xmin is the minimum and xmax is the maximum value of the property variable. We set x1/2 = xmin ,
xN+1/2 = xmax ,
xi+1/2 = xmin + i xi
0
(2.11)
for i = 1, 2, . . . , N − 1.
(3.1)
The product rule implies Furthermore we set
∂[G(t, x)f (t, x)] ∂[xG(t, x)f (t, x)] =x + G(t, x)f (t, x). ∂x ∂x (2.8)
xi =
(xi−1/2 + xi+1/2 ) , 2
xi = xi+1/2 − xi−1/2 .
(3.2)
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
3.1.2. Geometric grid In order to capture the initial profile properly, in some test problems it is useful to use a geometric grid. A typical onedimensional geometric grid is given as x1/2 = xmin ,
xi+1/2 = xmin + 2(i−N)/q (xmax − xmin )
for i = 1, 2, . . . , N
(3.3)
with xi and xi as given in (3.2) and the parameter q is any positive integer. Let Ωi = xi−1/2 , xi+1/2 for i ≥ 0. We approximate the initial data f˜ 0 (x) in each grid cell by ˜fi (0) = 1 f˜ 0 (x) dx. (3.4) xi Ωi 3.2. Numerical methods After having a discretized computational domain and assigning the initial data to each grid cell, the next step is to solve the given PBE (2.11) and (2.12). Here, we give two approaches for this purpose. 3.2.1. Method I: combination of FVS and MOC In this case we are interested to apply the method of characteristics for the growth process and a finite volume scheme for the aggregation process. For a scalar linear conservation law, for example PBE in the present case, there exist unique characteristic curves along which information propagates. If the solution moves along the propagation path-line, the advection term ∂Gf˜ /∂x in the current PBE disappears. 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 MOC we use (3.1) or (3.3) 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(t, x). dt Then Eq. (2.11) implies ∂ dx ˜ ∂f˜ (t, x) + f (t, x) ∂t ∂x dt
(3.5)
∂F(t, x) G(t, x)f˜ (t, x) ˜ =− + + S(t, x). (3.6) ∂x x Integration over the control volume Ωi (t) = [xi−1/2 (t), xi+1/2 (t)] implies xi+1/2 (t) dx ˜ ∂f˜ (t, x) dx + f (t, x) ∂t dt Ωi (t) xi−1/2 (t) G(t, x)f˜ (t, x) = − Fi+1/2 − Fi−1/2 + dx x Ωi (t) ˜ x) dx. S(t, (3.7) + Ωi (t)
1579
By using the Leibniz formula (Abramowitz & Stegun, 1972) backwards on the left hand side of (3.7), we get d f˜ (t, x) dx = − Fi+(1/2) − Fi−(1/2) dt Ωi (t) G(t, x)f˜ (t, x) ˜ x) dx. S(t, + dx+ x Ωi (t) Ωi (t) (3.8) Let f˜ i = f˜ i (t) and S˜ i = S˜ i (t) denote, respectively, the average values of the number density and nucleation term in each cell Ωi . They are defined as ˜ x) dx. ˜Si := 1 ˜fi := 1 ˜f (t, x) dx, S(t, xi (t) Ωi (t) xi (t) Ωi (t) (3.9) After using the above definitions, Eq. (3.8) implies d (xi+1/2 (t) − xi−1/2 (t))f˜ i dt Gi+1/2 f˜ i = − Fi+1/2 − Fi−1/2 + xi (t) + xi (t) S˜ i . xi (t) (3.10) By using the product rule and (3.5), the left hand side of (3.10) gives d (xi+1/2 (t) − xi−1/2 (t))f˜ i dt df˜ i dxi+1/2 dxi−1/2 ˜ = xi (t) + − fi dt dt dt = xi (t)
df˜ i + (Gi+1/2 − Gi−1/2 )f˜ i . dt
(3.11)
After replacing the left hand side of (3.10) with (3.11) and dividing the resultant equation by xi (t) one gets df˜ i 1 =− Fi+1/2 − Fi−1/2 dt xi (t) +
Gi+1/2 f˜ i f˜ i − (Gi+1/2 − Gi−1/2 ) + S˜ i , xi (t) xi (t)
(3.12)
where the cell interface values Fi+1/2 according to Filbet and Laurenc¸ot (2004) are given as ⎧ ⎪ i N ⎨ β(x , xk ) f˜ j h xk (t) f˜ k dx Fi+1/2 = ⎪ x ⎩j=αi,k Λ k=0 j
αi,k−1/2
+ f˜ αi,k −1 xi+1/2 −xk
β(x , xk ) x
dx
⎫ ⎪ ⎬ ⎪ ⎭
.
(3.13)
1580
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
Here, the integer αi,k corresponds to the index of the cell such that xi+1/2 (t) − xk (t) ∈ Ωαi,k −1 (t). In summary, we have to solve the following set of equations:
(3.14)
Integration of (2.11) over the control volume Ωi = [xi−1/2 , xi+1/2 ] implies ∂f˜ (t, x) ∂[G(t, x)f˜ (t, x)] G(t, x)f˜ (t, x) dx+ dx− dx ∂t ∂x x Ωi Ωi Ωi ∂F(t, x) =− S˜ i dx. (3.20) dx + ∂x Ωi Ωi
(3.15)
Let f˜ i = f˜ i (t) and S˜ i = S˜ i (t) denote respectively the average value of the number density and the nucleation term as given in (3.9). Then (3.20) becomes
df˜ i 1 =− Fi+1/2 − Fi−1/2 dt xi (t) +
Gi+1/2 f˜ i f˜ i − (Gi+1/2 − Gi−1/2 ) + S˜ i , xi (t) xi (t)
d˜xi+1/2 = Gi+1/2 dt
∀ i = 1, 2, . . . , N
with initial condition: f˜ (0, xi ) = f˜ 0 (xi ).
(3.16)
As a boundary condition the number density is taken to be zero outside the computational domain. This means that there exist no smaller or bigger particle than our specified size range of particles. The above system of ordinary differential equations can be solved by any standard ODE solver. 3.2.2. Presence of a nucleation term In order to overcome the nucleation problem 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 f˜ i (t) and xi (t) are initiated at these time levels and the time integrator restarts. In a special case where the nucleation only take place at the minimum particle size xmin , that is S(t, x) = B0 (t)δ(x − xmin ), one can use the nucleation simply as a boundary condition. Here, B0 (t) is the nucleation rate of particles of minimum particle size xmin and δ is a Dirac delta function, see Hounslow (1990). In that case (3.14)–(3.16) modifies to df˜ i 1 =− Fi+1/2 − Fi−1/2 dt x Gi+1/2 f˜ i f˜ i + − (Gi+1/2 − Gi−1/2 ) , xi xi d˜xi+1/2 = Gi+1/2 dt
∀ i = 1, 2, . . . , N
(3.17)
(3.18)
with initial and boundary conditions f˜ (0, xi ) = f˜ 0 (xi ),
xmin B0 (t) . f˜ (t, xmin ) = G(xmin )
(3.19)
3.2.3. Method II: semi-discrete finite volume scheme (FVS) Alternatively, one can use a semi-discrete finite volume scheme in order to solve the PBE (2.11) without any further modification. Here, we will give only the direct formulations of the scheme. For a detailed study of the scheme, the reader is referred to our articles (Qamar et al., 2006, 2007; Qamar & Warnecke, 2006).
∂f˜i 1 ˜ =− (Gf )i+1/2 − (Gf˜ )i−1/2 ∂t x Gi+1/2 f˜ i 1 − Fi+1/2 − Fi−1/2 + + S˜ i , x xi
(3.21)
where Fi+1/2 is given by (3.13). The flux (Gf˜ )i+1/2 at the right cell interface is given by the following flux-limited formula: xi (Gf˜ )i+1/2 = Gi+1/2 f˜ i + Φ(ri+ )(f˜ i+1 − f˜ i ) , 2 xi−1/2 (3.22) where xi−1/2 = xi − xi−1 ,
xi = xi+1/2 − xi−1/2 ,
(3.23)
and the flux limiting function Φ according to van Leer (1985) is defined as Φ(ri+ ) =
|ri+ | + ri+ . 1 + |ri+ |
(3.24)
The argument ri+ of the function Φ is the so-called upwind ratio of two consecutive solution gradients ri+ =
f˜ i − f˜ i−1 + ε . f˜ i+1 − f˜ i + ε
(3.25)
Analogously, one can formulate the flux (Gf˜ )i−1/2 at the left boundary of the control volume Ωi . The expression has to be evaluated with a small parameter, for example ε = 10−10 , in order to avoid division by zero. There are several other limiting functions proposed in the literature, namely, minmod, superbee and MC limiters, and so forth. Each of these limiters leads to a different high-resolution scheme, see LeVeque (2002). The above flux-limited formulation for the scheme guarantees the monotonicity and avoids the possible occurrence of wiggles and negative solution values. In case of discontinuous solutions the scheme gives second-order accuracy in smooth regions while first-order accuracy in the vicinity of shock discontinuities. This formula cannot be applied up to and including boundaries. Therefore, one has to use the first-order upwind formulation in the first two left boundary cells and in the last cell at the right boundary. Again, the resultant ODEs can be solved with any standard ODE-solver, for example the RK45 method in our case which is an embedded Runge-Kutta method of order four and five.
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
4. Numerical test problems The numerical schemes derived in this article can be applied to solve various combinations of nucleation, growth and aggregation processes for different choices of nucleation, growth and aggregation rates. However, to validate the current schemes, it is needed to test the schemes for the test problems where analytical solutions are also available. In the literature, analytical solutions of the PBEs are only available for the following combinations of processes: • Pure growth. • Simultaneous growth and aggregation. • Simultaneous growth and nucleation. In this section, we will consider test cases for above combinations of processes. Note that, analytical solutions for pure aggregation also do exist. However, this case was thoroughly discussed in Filbet and Laurenc¸ot (2004) and Qamar and Warnecke (2006) and will be omitted here. In all numerical test problems the lines represent the analytical solutions while symbols are used to represent the numerical solutions. 4.1. Pure growth 4.1.1. Constant growth problem 1 Here we consider constant growth with G(x) = G0 and an exponential initial distribution of the form: N0 x , (4.1) exp − f (0, x) := f0 (x) = x0 x0 where N0 and x0 are the total number of particles and the mean volume of particles at time t = 0, respectively. In both methods we use a geometric grid (3.3) with q = 3. The analytical solution is given by f (t, x) = f0 (x − G0 t).
(4.2)
Here, we take G0 = 1, N0 = 5, x0 = 0.01 and 60 mesh points. In Fig. 1, the numerical results of both methods are compared with the analytical solution. Note that the first plot which shows the MOC results uses a log scale on both axes so that the solutions are not merely shifted in the graph. The MOC results are in excellent agreement with the analytical solution and are hardly distinguishable. The stability and numerical diffusion of the numerical solution which marred the performance of the finite difference techniques are completely absent in the present scheme. The numerical results of FVS are also in agreement with the analytical solution. However, one can see a clear numerical dissipation in the solution. A comparison of the MOC and FVS shows that FVS results are little bit diffusive. Instead of any other scheme which uses some discretization technique for the advection term, the MOC results show a clear advantage where the advection term disappears from the governing equation. The first plot in Fig. 1 shows that the rapid growth rate results in very rapid growth of small particles, while the large particles stay nearly unchanged.
Fig. 1. A comparison of the results for constant growth problem 1.
1581
1582
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
Fig. 2. A comparison of the results for constant growth problem 2.
4.1.2. Constant growth problem 2 We consider a constant growth problem with initial data of the form: ⎧ 1 ⎪ ⎪ √ exp(−500(0.1x−0.3)2 ) if 2< x ≤4, ⎪ ⎪ ⎪ 0.32π ⎪ ⎪ ⎪ ⎨1 if 6< x ≤ 8, f (0, x)= 1−|x−11| if 10< x ≤ 12, ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 1−100(0.1x−1.5)2 if 14< x ≤ 16, ⎪ ⎪ ⎩ 0 elsewhere. We take the growth rate G(x) = 0.1 and divide the volume xmax = 20 into 100 equal subintervals. In Fig. 2, the numerical results of both methods are compared with the analytical solution and a first-order finite volume scheme. The MOC results are in excellent agreement with the analytical solution and are hardly distinguishable. The numerical results of FVS are also in agreement with the analytical solution. However, one can see a clear numerical dissipation in the solution. Finally, the firstorder finite volume scheme results are more diffusive than both schemes. 4.1.3. Linear growth Here G(x) = G0 x with exponential initial distribution (4.1). For the numerical simulation we take G0 = 1, x0 = 0.01 and N0 = 5. The analytical solution is given by Kumar and Ramkrishna (1997): f (t, x) = f0 (x e−G0 t ) e−G0 t .
(4.3)
In both methods we use a geometric grid (3.3) with q = 4 and 200 mesh points. Fig. 3 shows a comparison of numerical and analytical results. Once again the numerical results of MOC are in excellent agreement with the analytical results. However, one can see a visible under estimation in the results of the FVS on right ends of the number density profiles at t = 4 and 10. The numerical results of the finite volume scheme also show numerical dissipation on the left ends of the number density profiles. The numerical results of the first scheme again justifies the use of MOC for growth term.
Fig. 3. A comparison of analytical and numerical results for linear growth problem.
4.2. Simultaneous growth and aggregation In this case, five analytical solutions are available from Ramabhadran, Peterson, and Seinfeld (1976). The cases represent various combinations of constant and linear growth rates, constant and sum aggregation kernels, as well as exponential and gamma initial distributions. Kumar and Ramkrishna (1997) have also considered all these test problems. The specific choices for each case are given in Table 1 and their corresponding analytical solutions are given in Appendix A. Here we have tested our schemes for all five cases. In each case, we take N0 = 5, x0 = 0.01. A geometric grid with q = 6 is chosen to cover the size range of interest at initial time. Simulation results are presented for different simulation times. In both methods we use 100 mesh points for all five test cases. In the following we discuss all these case one by one. 4.2.1. Constant growth and constant aggregation Here we consider G(x) = G0 and β(x, x ) = β0 . The simulation results for cases 1 and 2 of Table 1 are shown in Figs. 4 and 5, respectively. Both figures indicate that the numerical results of
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
1583
Table 1 Various combinations of growth function, aggregation kernel and initial conditions tested for simultaneous growth and aggregation Case
G(x)
β(x, x )
1
1
100
2
1
100
3
x
10
4
x
10
5
x
x + x
f0 (x)
N0 x exp − x0 x0 x N0 x exp − x0 x0 x x 0 N0 exp − x0 x 0 x N0 x exp − x0 x0 x 0 x N0 exp − x0 x0
the first and second methods are in very good agreement with the analytical results and those in Kumar and Ramkrishna (1997) across several order of magnitude. Even at very small values the results looks very promising. However, the numerical results of
Fig. 5. A comparison of analytical and numerical results for simultaneous growth and aggregation, case 2 (Table 1).
Fig. 4. A comparison of analytical and numerical results for simultaneous growth and aggregation, case 1 (Table 1).
the second method are dissipative at the left discontinuity of the number density profile and is in good agreement in the other regions. There are some differences between first scheme results and analytical and numerical results on the far left side of the number density profile for t = 0.01 in Fig. 5 and similar but not visible differences in all other simulations in Fig. 4 are because the analytical solution in this range breaks down, see also Kumar and Ramkrishna (1997). Ramabhadran et al. (1976) have pointed out that for approximately equal growth and aggregation rates, their analytical solution for G(x) = G0 and β(x, x ) = β0 is valid for only large particle sizes. For the aggregation-dominated situation, the range of validity is increased, where as for the growth-dominated situation, the solution breaks down completely. The simulations given in Figs. 4 and 5 are somehow aggregation-dominated cases in order to allow analytical solutions to hold in a reasonably wide range. For this problem
1584
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
Fig. 7. A comparison of analytical and numerical results for simultaneous growth and aggregation, case 4 (Table 1).
Fig. 6. A comparison of analytical and numerical results for simultaneous growth and aggregation, case 3 (Table 1).
with exponential initial distribution the extent of aggregation M0 (t)/M0 (0), as growth processes conserve total number, is 0.288 and 3.998 × 10−4 at times t = 10−2 and 10, respectively. Similarly the extent of growth processes M1 (t)/M1 (0), as the aggregation process conserve mass, is 1.501 and 4.13 at times t = 10−2 and 10, respectively. 4.2.2. Linear growth and constant aggregation In this case we consider G(x) = G0 x and β(x, x ) = β0 . Figs. 6 and 7 show evolving size distributions for cases 3 and 4 of Table 1. The numerical results are once again in very good agreement with the analytical results and those in Kumar and Ramkrishna (1997). Again, there is a visible smearing in the second method results at far left end of the number density profiles, while are in good agreement with analytical and MOC solutions
in other regions. In case of exponential initial distribution the extent of aggregation M0 (t)/M0 (0) is 0.039 and 3.98 × 10−3 at times t = 1 and 10, respectively. Similarly the extent of growth processes M1 (t)/M1 (0) is 2.72 and 2.2 × 104 at times t = 1 and 10, respectively. 4.2.3. Linear growth and sum aggregation Here we consider G(x) = G0 x and β(x, x ) = β0 (x + x ). The numerical results for the case 5 of Table 1 are presented in Fig. 8. It shows that even with size dependent rate functions, the analytical and the numerical results of first and second methods are in excellent agreement. However, in the results of the second method there is a visible smearing at the left end of the number density profiles and very slight over prediction at the right end. The above test simulations clearly show that the proposed techniques are quite robust and yield very good results for simultaneous growth and aggregation. Here the extent of aggregation M0 (t)/M0 (0) is 0.92 and 0.39 at times t = 1 and 3, respectively. Similarly extent of growth processes
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
M1 (t)/M1 (0) is 2.72 and 20.09 at times t = 1 and 3, respectively. 4.3. Simultaneous nucleation and growth In this section we consider some numerical test problems in order to test the proposed schemes for simultaneous nucleation and growth processes. Most of these numerical problems were also considered by Kumar and Ramkrishna (1997) in order to test their scheme. In both methods we use again geometric grid with q = 6. We take 100 mesh points in both numerical examples.
The volume and time ranges are given by 0 ≤ x ≤ 2.0 and 0 ≤ t ≤ 0.5, respectively. The square step initial condition for the number density is given as (Lim et al., 2002): 100 if 0.4 < x ≤ 0.6, f (0, x) = 0.01 elsewhere. Here we consider a constant growth with G = 1.0 and a linear grid of 200 mesh points. The analytical solution is given as f (t, x) =
4.3.1. Constant growth stiff nucleation This problem was considered by Lim et al. (2002). Suppose that the stiff nucleation takes place at a minimum volume as a function of time: f (0, t) = 100 + 106 exp(−104 (t − 0.215)2 ).
1585
⎧ 2 ⎨ 100 − 106 exp(−104 (Gt − x − 0.215) ) if 0.0 < x ≤ Gt, ⎩
100
if 0.4 < x − Gt ≤ 0.6,
0.01
elsewhere.
The results are shown in Fig. 9. The numerical results of both methods are compared with the analytical solutions and firstorder finite volume scheme. The MOC results are in excellent agreement with the analytical solution and are hardly distinguishable. The numerical results of FVS are diffusive as compared to MOC. Finally, the first-order finite volume scheme results are more diffusive than both schemes. 4.3.2. Constant growth and exponential nucleation In this case we take G(x) = G0 , and S(x) = B0 /x0,n exp(−x/x0,n ). As an initial condition we take the exponential initial distribution of the form f0 (x) = N0 /x0 exp(−x/x0 ). For numerical simulation we take G0 = 1, N0 = 10, x0 = 0.01, x0,n = 0.001, and B0 = 105 . The analytical solution for this case is given as (Kumar and Ramkrishna, 1997): xlow x B0 exp − − exp − , f (t, x)=f0 (x − Gt) + G0 x0,n x0,n (4.4)
Fig. 8. A comparison of analytical and numerical results for simultaneous growth and aggregation, case 5 (Table 1).
Fig. 9. A comparison of exact and numerical results for simultaneous nucleation and growth.
1586
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
Fig. 11. A comparison of the results for linear growth and exponential nucleation. Fig. 10. A comparison of the results for constant growth and exponential nucleation.
and Ramkrishna, 1997): f (t, x) = f0 (x e−G0 t ) exp(−G0 t) +
xlow = max(x1 , x − G0 t).
(4.5)
Here x1 represent the center of the first cell in the mesh. Since now nucleation takes place, we have to add a new cell of nuclei size at each time step. In order to keep the number of cells fixed we delete the end cell at the same time level. The numerical results are shown in Fig. 10. Again the numerical results of both methods are in very good agreement with the analytical solutions. The results of the second method are also excellent with a minor dissipation in the numerical result at t = 10−2 .
− exp −
x x0,n
,
xlow B0 exp − G0 x x0,n
xlow = max (x1 , x e−G0 t ) . (4.6)
The numerical results are shown in Fig. 11. The results of the MOC are in excellent agreement with the analytical solution. However, there is visible dissipation in the numerical results of the FVS. From this and previous results one can conclude that MOC become more important when the growth term is not constant. 4.4. Simultaneous nucleation, growth and aggregation
4.3.3. Linear growth and exponential nucleation Here the initial data and the nucleation term have the same formulation as given in the last example but now G(x) = G0 x. Here we take G0 = 1, N0 = 10, x0 = 0.01, x0,n = 0.001, and B0 = 10. The analytical solution for this case is given as (Kumar
Unfortunately, there is no analytical solution available in this case. We consider the constant growth rate G(x) = G0 , exponential nucleation rate S(x) = B0 /x0,n exp(−x/x0,n ) and a constant aggregation kernel β(x, x ) = β0 . As an initial con-
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
1587
5. Conclusions
Fig. 12. A comparison of FVS and MOC + FVS results with each other and analytical results for nucleation, growth and agglomeration.
dition we take the exponential initial distribution of the form f0 (x) = N0 /x0 exp(−x/x0 ). For numerical simulation we take G0 = 1, N0 = 10, x0 = 0.01, B0 = 1, x0,n = 0.001, β0 = 100, and 100 mesh points. The numerical results for both methods are compared with each other in Fig. 12. The second plot in Fig. 12 represents the results of both methods in the absence of aggregation. The first plot shows that first method resolve the discontinuities very well as compared to the second method. A comparison of both plots, that is with and without aggregation, shows that aggregation significantly effects the distribution of particles present initially and those born afterwards. Particularly, the distribution of the nucleated particles, which lie on the left hand side of discontinuity in the size range, is changed completely in the presence of aggregation. In the second plot since the analytical solutions and MOC solutions were completely overlapping, therefore we omit the analytical solution here. The numerical results of both methods are in good agreement with each other. However, the results of the first method are more resolved as compared to the second method where dissipation is visible at the steep discontinuities.
In this article two numerical methods are proposed for the solution of PBEs for simultaneous nucleation, growth and aggregation processes. The first method combines the finite volume scheme (FVS) for aggregation process with a method of characteristics (MOC) for growth process. For handling nucleation term in this method, a cell of nuclei size is added at a given time level. The second method purely uses a finite volume scheme for nucleation, growth and aggregation of particles. Both methods use the same finite volume scheme for the aggregation process. From the numerical test problems, it was found that the use of MOC for linear growth processes is more successful as compared to any finite difference and even finite volume schemes. Although the growth process looks simpler than the aggregation process, still it give rise to complexity in the PBE. Since the presence of growth term introduces hyperbolicity in the governing equation, there may exist discontinuities or steep fronts in the solution. The approximation of advection term for the growth process by a finite difference or a finite volume scheme introduces an inherent viscosity in the solution. This leads to the smearing of the resultant solution in the vicinity of discontinuities (steep fronts). In case of MOC, the solution moves along the propagation path-line and the advection term ∂Gf˜ /∂x simply disappears from the governing equation. Therefore, no special care is needed for the advection term which is usually required in other schemes. Even though the results of MOC are superior to the FVS, still the FVS worked very well for all numerical test cases presented in this article. The use of FVS for the current integrodifferential population balance equation show their applicability, generality and robustness. The schemes also show their flexibility to be adapted in the same manner and to the same grids which were basically proposed for the methods specific to PBEs. Note that, the aggregation term in both schemes is calculated with the same finite volume scheme. This shows that finite volume scheme works very well for aggregation terms and gives very accurate results as for as the advection terms is approximated very well. Since, in the first method the growth term is calculated with the MOC, which gives more accurate results, the overall results for simultaneous growth and aggregation are also very excellent and almost indistinguishable from the analytical results. Especially, for the two-component problems involving simultaneous growth and aggregation processes, the use of pure finite volume seems to be computationally expensive. It was found that, the use of MOC for growth terms makes the first method at least two times faster than the second method. However, in most of the test cases presented in this article the computational time was not more than a minute for both methods. There were only one or two examples where the computational time was a little bit more than a minute. Several numerical results in different combinations of nucleation, growth and aggregation are presented in this article. The numerical results of the schemes were found to be in good agreement with the analytical solutions.
1588
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
Here Λ = G0 /β0 N0 x0 and Q + Q2 − R g1 = − , Rx0
Table 2 Moments M0 and M1 of the distributions Case
Moments
β(x, x ) = β0
M0 =
G(x) = G0
M1 = N0 x0 1 −
β(x, x ) = β0
M0 =
G(x) = G0 x β(x, x )
= β0 (x
2N0 2 + β0 N0 t
2G0 ln β0 N0 x0
2 2 + β0 N0 t
2N0 2 + β0 N0 t
M1 = N0 x0 exp(β0 t) + x )
G(x) = G0 x
M0 = N0 exp
β N x 0 0 0
G0 M1 = N0 x0 exp(β0 t)
(1 − exp(G0 t))
g2 = −
Q(x) = χ − Λ(χ − 1 − χ ln χ),
Q−
χ=
Q2 − R , Rx0
N0 , M0
R(χ) = χ − (2Λ)2 χ(χ − 1 − ln χ) − (χ − 1)2 .
(A.5)
A.2. Linear growth and constant aggregation In case of linear growth G(x) = G0 x and size-independent aggregation kernel β(x, x ) = β0 , the analytical solution for exponential initial distribution (A.1) is given by M02 M0 f (t, x) = exp − x . (A.6) M1 M1
Appendix A. Analytical solutions for simultaneous growth and aggregation problems There are several combinations of growth rates and aggregation kernels where analytical solutions are available. These solution were derived by Ramabhadran et al. (1976). Since two simultaneous processes are taking place, the dynamics of the particle size distribution depends on two characteristic times, one for aggregation and other for growth. The ratio of these two times is the basic dimensionless parameter of the problem. Let us define Λ as the ratio of the characteristic times of growth and aggregation. The analytical solution for two initial distribution and different combinations of the growth and aggregation rates are given below.
For Gaussian-type initial distribution (A.3) the analytical solution is given as f (t, x) =
M02 1 √ M1 1 − M0 /N0 N0 x M 0 N0 x . sinh × exp − 1− M1 N0 M1
(A.7)
The moments M0 and M1 for this case are given in Table 2. A.3. Linear growth and sum aggregation
A.1. Constant growth and constant aggregation For the constant growth G(x) = G0 and size-independent aggregation kernel β(x, x ) = β0 , the analytical solution with exponential initial distribution: N0 x , (A.1) exp − f (0, x) := f0 (x) = x0 x0 is given as f (t, x) =
M02 /M1 1 − 2Λx0 (N0 − M0 /M1 ) M0 (x − 2Λx0 (N0 /M0 − 1)) × exp − , M1 1 − 2Λx0 (N0 − M0 /M1 )
For the case that G(x) = G0 x and β(x, x ) = β0 (x + x ), the analytical solution for exponential initial distribution (A.1) is given as f (t, x) =
M0 √ 1 − M0 /N0 M0 2N0 M0 N0 x , × exp − −1 x I1 2 1− M1 M0 N0 M1 x
(A.8) (A.2)
where Λ = G0 /β0 N0 x0 . The moments M0 and M1 are given in Table 2. For Gaussian-type initial distribution N0 x x , (A.3) f (0, x) := f0 (x) = 2 exp − x0 x0 the analytical solution is given by M0 /x02 N0 f (t, x) = −1 exp g1 x − 2Λx0 R(g1 − g2 ) M0 N0 − exp g2 x − 2Λx0 −1 . (A.4) M0
where I1 is the modified Bessel function of first kind of order one. The moments M0 and M1 for are given in Table 2. References Abramowitz, M., & Stegun, I.A. (1972). Handbook of mathematical functions with formulas, graphs and mathematical tables (p. 11). New York: Dover. David, R., Villermaux, P., & Klein, J. P. (1991). Crystallization and precipitation engineering-IV. Kinetic model to adipic acid crystallization. Chemical Engineering Science, 46, 1129–1136. Filbet, F., & Laurenc¸ot, P. (2004). Numerical simulation of the Smoluchowski coagulation equation. SIAM Journal on Scientific Computing, 25(6), 2004–2048. Gunawan, R., Fusman, I., & Braatz, R. D. (2004). High resolution algorithms for multidimensional population balance equations. AIChE Journal, 50, 2738–2749. Hounslow, M. J. (1990). A discretized population balance for continuous systems at steady-state. AIChE Journal, 36, 106–116.
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589 Hounslow, M. J., Ryall, R. L., & Marshall, V. R. (1988). A discretized population balance for nucleation, growth, and aggregation. AIChE Journal, 34, 1821–1832. 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. Lapidus, L., & Pinder, G. F. (1982). Numerical solutions of partial differential equations in science and engineering. New York, USA: Wiley. LeVeque, R. J. (2002). Finite volume methods for hyperbolic problems. New York, NY: Cambridge University Press. 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. Litster, J. D., Smith, D. J., & Hounslow, M. J. (1995). Adjustable discretized population balance for growth and agglomeration. AIChE Journal, 41, 591–603. 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 and Chemical Engineering, 26, 1103–1116. Motz, S., Mitrovi´c, A., & Gilles, E. D. (2002). Comparison of numerical methods for the simulation of dispersed phase systems. Chemical Engineering Science, 57, 4329–4344. Muhr, H., David, J., Villermaux, J., & Jezequel, P. H. (1996). Crystallization and precipitation engineering VI: Solving population balance in the case of the precipitation of silver bromide crystals with high primary nucleation rate by
1589
using first order upwind differentiation. Chemical Engineering Science, 51, 309–319. Qamar, S., Ashfaq, A., Warnecke, G., Angelov, I., Elsner, M. P., & SeidelMorgenstern, A. (2007). Adaptive high resolution schemes for population balances in crystallization processes. Computers and Chemical Engineering, 31, 1296–1311. 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 and Chemical Engineering, 30, 1119–1131. Qamar, S., & Warnecke, G. (2006). Solving population balance equations for two-component aggregation by a finite volume scheme. Chemical Engineering Science, 62, 679–693. Ramabhadran, T. E., Peterson, T. W., & Seinfeld, J. H. (1976). Dynamics of aerosol coagulation and condensation. AIChE Journal, 22, 840– 851. Ramkrishna, D. (1985). The status of population balances. Reviews in Chemical Engineering, 3(1), 49–95. 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. 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. van Leer, B. (1985). Upwind-difference methods for aerodynamic problems governed by the Euler equations. Lectures in applied mathematics (Vol. 22-Part 2, pp. 327–336). Providence, RI: American Mathematical Society.
Numerical solution of population balance equations for nucleation, growth and aggregation processes Shamsul Qamar ∗ , Gerald Warnecke Institute for Analysis and Numerics, Otto-von-Guericke University, PSF 4120, 39106 Magdeburg, Germany Received 26 August 2006; received in revised form 16 January 2007; accepted 16 January 2007 Available online 28 January 2007
Abstract This article focuses on the derivation of numerical schemes for solving population balance models (PBMs) with simultaneous nucleation, growth and aggregation processes. Two numerical methods are proposed for this purpose. The first method combines a method of characteristics (MOC) for growth process with a finite volume scheme (FVS) for aggregation process. For handling nucleation terms, a cell of nuclei size is added at a given time level. The second method purely uses a semi-discrete finite volume scheme for nucleation, growth and aggregation of particles. Note that both schemes use the same finite volume scheme for aggregation process. On one hand, the method of characteristics offers a technique which is in general a powerful tool for solving linear growth processes, has the capability to overcome numerical diffusion and dispersion, is computationally efficient, as well as give highly resolved solutions. On the other hand, the finite volume schemes which were derived for a general system in divergence form, are applicable to any grid to control resolution, and are also computationally not expensive. In the first method a combination of finite volume scheme and the method of characteristics gives a highly accurate and efficient scheme for simultaneous nucleation, growth and aggregation processes. The second method demonstrates the applicability, generality, robustness and efficiency of high-resolution schemes. The proposed techniques are tested for pure growth, simultaneous growth and aggregation, nucleation and growth, as well as simultaneous nucleation, growth and aggregation processes. The numerical results of both schemes are compared with each other and are also validated against available analytical solutions. The numerical results of the schemes are in good agreement with the analytical solutions. © 2007 Elsevier Ltd. All rights reserved. Keywords: Population balance models; High-resolution finite volume schemes; Method of characteristics; Nucleation; Growth; Aggregation; Hyperbolic conservation laws; Particles
1. Introduction The dynamic behavior of a population of small particles is a subject of interest in various engineering fields. These areas include atmospheric physics, precipitation, crystallization, pharmaceutical manufacture, aerosol formation, colloid chemistry, growth of microbial and cell populations, and so on. Most of such systems involve simultaneous nucleation, growth and aggregation of particles. Population balances, in their integrodifferential equation form, are a widely used tool for simulating these processes. Since, population balance equations (PBEs) can only be solved exactly in simplified cases, numerical solutions are usually needed. Many studies are therefore focused on
∗
Corresponding author. Tel.: +49 391 6712027; fax: +49 391 6718073. E-mail address: [email protected] (S. Qamar). 0098-1354/$ – see front matter © 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2007.01.006
the accurate and efficient numerical solutions of the population balance equations. In the past several numerical techniques were introduced for solving PBEs, see for example Ramkrishna (1985) for a review. Discussion on the stability and numerical dispersion of various finite difference approximations can be found in Lapidus and Pinder (1982). An extension of the finite difference-type discretization methods have been proposed by Hounslow, Ryall, and Marshall (1998) and David, Villermaux, and Klein (1991). However, their numerical results were found to be unsatisfactory, see Kumar and Ramkrishna (1997). A second-order spatial discretization for the growth term as well as an accurate and complex discretization for the aggregation term on a geometric grid was proposed by Litster, Smith, and Hounslow (1994). Later on, Muhr, David, Villermaux, and Jezequel (1996) used a firstorder upwind scheme for the spatial discretization of the growth term. However, due to numerical diffusion a loss in accuracy near steep fronts or discontinuities was observed.
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
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 term in order to solve PBEs with simultaneous nucleation, growth and aggregation processes. They found that, in contrast to other numerical techniques, MOC avoids the numerical dissipation error caused by the growth term discretization. For handling the nucleation term, which is usually difficult to treat, MOC was combined with a procedure of adding a cell of the nuclei size at a given time level. A standard ordinary differential equation (ODE) solver was used to solve the resultant ODEs. A short overview of the methods of Hounslow et al. (1988) and David et al. (1991) is also given in their article. The finite volume schemes, which were originally developed for the gas dynamics, were also applied for the numerical solution of population balance models in crystallization process, see Gunawan, Fusman, and Braatz (2004) and Ma, Tafti, and Braatz (2002a, 2002b) and Motz, Mitrovi´c, and Gilles (2002). Further, Lim, Lann, and Meyer (2002) applied highresolution spatial discretization methods (WENO schemes) and the modified method of characteristics (MOC) for the dynamic simulations of the batch crystallization including nucleation, growth, aggregation and breakage kinetics. They observed that, steep moving fronts or discontinuities appearing in the solution of the crystallization models were captured well by both methods without any spurious oscillations. Apart from these methods, Filbet and Laurenc¸ot (2004) proposed a numerical scheme for single-component aggregation which is based on a conservative finite volume formulation. The authors have rewritten the population balance equation for aggregation problems in a form which can be readily solved by a finite volume scheme. This special reformulation was a great achievement which enables one to apply the finite volume scheme in the aggregation case as well. Later on, Qamar and Warnecke (2006) have extended the scheme in order to simulate two-component aggregation problems. Apart from that, Qamar, Elsner, and Angelov (2006) and Qamar et al. (2007) have applied high-resolution finite volume schemes for solving population balance models in crystallization processes, where they considered growth and nucleation kinetics while neglecting the aggregation and breakage processes. The aim of the current work is to numerically solve PBEs for simultaneous nucleation, growth and aggregation processes. Two methods are proposed in this article. The first method combines a method of characteristics for growth term (Kumar & Ramkrishna, 1997) and a finite volume scheme (Filbet & Laurenc¸ot, 2004; Qamar & Ashfaq, 2006) for aggregation term. In the second method, a semi-discrete finite volume scheme (FVS) is used for all processes, see Qamar et al. (2006, 2007) and Qamar and Warnecke (2006) for a review. Note that, in both methods the aggregation kinetics is solved with the same finite volume formulation. The main difference in both methods is the way they calculate the growth term, for example the first method uses MOC for the growth term and the advection term simply disappears from the governing equation, while the second method uses a finite volume scheme to discretize the advection term. The efficiency and accuracy of the resultant methods are analyzed by comparing their numerical results with each other and with avail-
1577
able analytical solutions. For handling the nucleation term in the first scheme, MOC is combined with a procedure of adding a cell of the nuclei size at each time level. The same procedure was also used in Kumar and Ramkrishna (1997). A standard ODEs solver can be used to solve the resultant ODEs. There are two main differences between the first method and the one used in Kumar and Ramkrishna (1997). In our current method the aggregation term is treated with a finite volume scheme, while in Kumar and Ramkrishna (1997) the authors have used the fixed pivot technique for this purpose. Secondly, our numerical method uses a reformulated PBE instead of the original one. This reformulation is necessary in order to apply the finite volume scheme. In case of reformulated PBE, our numerical method calculates the volume (mass) density instead of the number density. However, one can easily recover back the number density at the end of simulation. Similarly Lim et al. (2002) used the original PBE and a WENO scheme which is far different from the current finite volume schemes. Note 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 method of characteristics to deal with the hyperbolic nature of the growth term, are needed when accurate population balance model solutions are needed to compute real-time 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 particle size distribution can be found in the recent articles of Shi, El-Farra, and Mhaskar (2006) and Shi, Mhaskar, and El-Farra (2005). Although we are not considering control problems in this work, the current schemes can be used in this direction as well. The paper is organized as follows. In Section 2, we start with an integro-differential PBE and reformulate it to a form which can be readily used for the derivation of our proposed numerical methods. In Section 3, we give two different grid discretizations and then derive both numerical methods for the solution of PBEs for simultaneous nucleation, growth and aggregation processes. In Section 4, we give some numerical results for different combinations of the nucleation, growth and aggregation processes. For validation of our numerical results we use the available analytical solutions. In Section 5, we give conclusions and remarks. 2. Mathematical model In this section we will reformulate the underlying onedimensional PBE in a form which can be readily used to derive the proposed numerical methods. Consider a PBE of the form: ∂f (t, x) ∂[G(t, x)f (t, x)] + ∂t ∂x = Qagg (t, x) + S(t, x),
(t, x) ∈ R2+ ,
(2.1)
1578
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
f (0, x) = f0 (x),
x ∈ R+ ,
(2.2)
where R+ :=]0, +∞[, f (t, x) is the number density; G(t, x) the growth term; S(t, x) the rate at which particle of size x are nucleated; Qagg (t, x) is the aggregation reaction term which is defined as x 1 Qagg (t, x) = β(t, x − x , x ) f (t, x − x ) f (t, x ) dx 2 0
∞ −
β(t, x, x ) f (t, x) f (t, x ) dx .
xi f (t, x) dx,
i = 0, 1, 2, . . . ,
(2.4)
0
where the zeroth moment M0 (t) and the first moment M1 (t) represent the total number of particles and the total volume of particles, respectively. We want to consider the mass balance. Multiplying both sides of Eqs. (2.1)–(2.3) by x, we get x
∂f (t, x) ∂[G(t, x)f (t, x)] +x ∂t ∂x = xQagg (t, x) + xS(t, x), x ∈ R+ ,
xf (0, x) = xf0 (x),
xQagg (t, x) =
1 2
x
(t, x) ∈ R2+ ,
(2.5) (2.6)
xβ(t, x − x , x ) f (t, x − x ) f (t, x ) dx
0
∞ −
xβ(t, x, x ) f (t, x) f (t, x ) dx .
(2.9)
where x ∞ u β(t, u, v) f (t, u) f (t, v) dv du,
F(t, x) = −
x ∈ R+ . (2.10)
˜ x) := xS(t, x). After Let us define f˜ (t, x) := xf (t, x) and S(t, using Eqs. (2.4)–(2.10) in (2.5) and (2.6), we obtain the mass balance: ∂f˜ (t, x) ∂[G(t, x)f˜ (t, x)] G(t, x)f˜ (t, x) + − ∂t ∂x x ∂F(t, x) ˜ =− + S(t, x), ∂x
f˜ (0, x) = f˜ 0 ,
(x, x ) ∈ R2+ .
It is basically the probability that particles of sizes x, x meet and stick together. The moments of the distribution function f are defined as Mi (t) =
∂F(t, x) , ∂x
0 x−u
The first term in Eq. (2.5) accounts for the formation of particles with volume x resulting from the merging of two particles with respective volume x and x − x , where x ∈ ]0, x[. The second term is the death term due to loss of particle of volumes x by aggregation with other particles of any size. The aggregation coefficient β = β(t, x, x ) characterizes the rate at which the aggregation of two particles with respective volumes x and x produces a particle of volume x + x and is a non-negative symmetric function,
∞
xQagg (t, x) = −
(2.3)
0
0 ≤ β(t, x, x ) = β(t, x , x),
Furthermore, the aggregation term xQagg (t, x) can be rewritten as (Filbet & Laurenc¸ot, 2004; Qamar & Ashfaq, 2006):
(2.7)
x ∈ R+ .
(2.12)
Instead of the original equation (2.1) we will discretize (2.11). Due to the relation f˜ (t, x) := xf (t, x) one can easily recover f (t, x). 3. Numerical procedure In the following we will derive two methods in order to solve the above modified population balance equations (2.11) and (2.12) for the simulation of simultaneous nucleation growth and aggregation processes. 3.1. Domain discretization In order to apply any numerical scheme, the first step is to discretize the computational domain. Here, we give two types of grid discretization, that is, regular (or irregular grid) and geometric grid. 3.1.1. Regular/irregular grid Let N be a large integer representing the total number of cells, we denote by (xi−(1/2) )i ∈ {1,...,N+1} a mesh of [xmin , xmax ], where xmin is the minimum and xmax is the maximum value of the property variable. We set x1/2 = xmin ,
xN+1/2 = xmax ,
xi+1/2 = xmin + i xi
0
(2.11)
for i = 1, 2, . . . , N − 1.
(3.1)
The product rule implies Furthermore we set
∂[G(t, x)f (t, x)] ∂[xG(t, x)f (t, x)] =x + G(t, x)f (t, x). ∂x ∂x (2.8)
xi =
(xi−1/2 + xi+1/2 ) , 2
xi = xi+1/2 − xi−1/2 .
(3.2)
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
3.1.2. Geometric grid In order to capture the initial profile properly, in some test problems it is useful to use a geometric grid. A typical onedimensional geometric grid is given as x1/2 = xmin ,
xi+1/2 = xmin + 2(i−N)/q (xmax − xmin )
for i = 1, 2, . . . , N
(3.3)
with xi and xi as given in (3.2) and the parameter q is any positive integer. Let Ωi = xi−1/2 , xi+1/2 for i ≥ 0. We approximate the initial data f˜ 0 (x) in each grid cell by ˜fi (0) = 1 f˜ 0 (x) dx. (3.4) xi Ωi 3.2. Numerical methods After having a discretized computational domain and assigning the initial data to each grid cell, the next step is to solve the given PBE (2.11) and (2.12). Here, we give two approaches for this purpose. 3.2.1. Method I: combination of FVS and MOC In this case we are interested to apply the method of characteristics for the growth process and a finite volume scheme for the aggregation process. For a scalar linear conservation law, for example PBE in the present case, there exist unique characteristic curves along which information propagates. If the solution moves along the propagation path-line, the advection term ∂Gf˜ /∂x in the current PBE disappears. 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 MOC we use (3.1) or (3.3) 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(t, x). dt Then Eq. (2.11) implies ∂ dx ˜ ∂f˜ (t, x) + f (t, x) ∂t ∂x dt
(3.5)
∂F(t, x) G(t, x)f˜ (t, x) ˜ =− + + S(t, x). (3.6) ∂x x Integration over the control volume Ωi (t) = [xi−1/2 (t), xi+1/2 (t)] implies xi+1/2 (t) dx ˜ ∂f˜ (t, x) dx + f (t, x) ∂t dt Ωi (t) xi−1/2 (t) G(t, x)f˜ (t, x) = − Fi+1/2 − Fi−1/2 + dx x Ωi (t) ˜ x) dx. S(t, (3.7) + Ωi (t)
1579
By using the Leibniz formula (Abramowitz & Stegun, 1972) backwards on the left hand side of (3.7), we get d f˜ (t, x) dx = − Fi+(1/2) − Fi−(1/2) dt Ωi (t) G(t, x)f˜ (t, x) ˜ x) dx. S(t, + dx+ x Ωi (t) Ωi (t) (3.8) Let f˜ i = f˜ i (t) and S˜ i = S˜ i (t) denote, respectively, the average values of the number density and nucleation term in each cell Ωi . They are defined as ˜ x) dx. ˜Si := 1 ˜fi := 1 ˜f (t, x) dx, S(t, xi (t) Ωi (t) xi (t) Ωi (t) (3.9) After using the above definitions, Eq. (3.8) implies d (xi+1/2 (t) − xi−1/2 (t))f˜ i dt Gi+1/2 f˜ i = − Fi+1/2 − Fi−1/2 + xi (t) + xi (t) S˜ i . xi (t) (3.10) By using the product rule and (3.5), the left hand side of (3.10) gives d (xi+1/2 (t) − xi−1/2 (t))f˜ i dt df˜ i dxi+1/2 dxi−1/2 ˜ = xi (t) + − fi dt dt dt = xi (t)
df˜ i + (Gi+1/2 − Gi−1/2 )f˜ i . dt
(3.11)
After replacing the left hand side of (3.10) with (3.11) and dividing the resultant equation by xi (t) one gets df˜ i 1 =− Fi+1/2 − Fi−1/2 dt xi (t) +
Gi+1/2 f˜ i f˜ i − (Gi+1/2 − Gi−1/2 ) + S˜ i , xi (t) xi (t)
(3.12)
where the cell interface values Fi+1/2 according to Filbet and Laurenc¸ot (2004) are given as ⎧ ⎪ i N ⎨ β(x , xk ) f˜ j h xk (t) f˜ k dx Fi+1/2 = ⎪ x ⎩j=αi,k Λ k=0 j
αi,k−1/2
+ f˜ αi,k −1 xi+1/2 −xk
β(x , xk ) x
dx
⎫ ⎪ ⎬ ⎪ ⎭
.
(3.13)
1580
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
Here, the integer αi,k corresponds to the index of the cell such that xi+1/2 (t) − xk (t) ∈ Ωαi,k −1 (t). In summary, we have to solve the following set of equations:
(3.14)
Integration of (2.11) over the control volume Ωi = [xi−1/2 , xi+1/2 ] implies ∂f˜ (t, x) ∂[G(t, x)f˜ (t, x)] G(t, x)f˜ (t, x) dx+ dx− dx ∂t ∂x x Ωi Ωi Ωi ∂F(t, x) =− S˜ i dx. (3.20) dx + ∂x Ωi Ωi
(3.15)
Let f˜ i = f˜ i (t) and S˜ i = S˜ i (t) denote respectively the average value of the number density and the nucleation term as given in (3.9). Then (3.20) becomes
df˜ i 1 =− Fi+1/2 − Fi−1/2 dt xi (t) +
Gi+1/2 f˜ i f˜ i − (Gi+1/2 − Gi−1/2 ) + S˜ i , xi (t) xi (t)
d˜xi+1/2 = Gi+1/2 dt
∀ i = 1, 2, . . . , N
with initial condition: f˜ (0, xi ) = f˜ 0 (xi ).
(3.16)
As a boundary condition the number density is taken to be zero outside the computational domain. This means that there exist no smaller or bigger particle than our specified size range of particles. The above system of ordinary differential equations can be solved by any standard ODE solver. 3.2.2. Presence of a nucleation term In order to overcome the nucleation problem 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 f˜ i (t) and xi (t) are initiated at these time levels and the time integrator restarts. In a special case where the nucleation only take place at the minimum particle size xmin , that is S(t, x) = B0 (t)δ(x − xmin ), one can use the nucleation simply as a boundary condition. Here, B0 (t) is the nucleation rate of particles of minimum particle size xmin and δ is a Dirac delta function, see Hounslow (1990). In that case (3.14)–(3.16) modifies to df˜ i 1 =− Fi+1/2 − Fi−1/2 dt x Gi+1/2 f˜ i f˜ i + − (Gi+1/2 − Gi−1/2 ) , xi xi d˜xi+1/2 = Gi+1/2 dt
∀ i = 1, 2, . . . , N
(3.17)
(3.18)
with initial and boundary conditions f˜ (0, xi ) = f˜ 0 (xi ),
xmin B0 (t) . f˜ (t, xmin ) = G(xmin )
(3.19)
3.2.3. Method II: semi-discrete finite volume scheme (FVS) Alternatively, one can use a semi-discrete finite volume scheme in order to solve the PBE (2.11) without any further modification. Here, we will give only the direct formulations of the scheme. For a detailed study of the scheme, the reader is referred to our articles (Qamar et al., 2006, 2007; Qamar & Warnecke, 2006).
∂f˜i 1 ˜ =− (Gf )i+1/2 − (Gf˜ )i−1/2 ∂t x Gi+1/2 f˜ i 1 − Fi+1/2 − Fi−1/2 + + S˜ i , x xi
(3.21)
where Fi+1/2 is given by (3.13). The flux (Gf˜ )i+1/2 at the right cell interface is given by the following flux-limited formula: xi (Gf˜ )i+1/2 = Gi+1/2 f˜ i + Φ(ri+ )(f˜ i+1 − f˜ i ) , 2 xi−1/2 (3.22) where xi−1/2 = xi − xi−1 ,
xi = xi+1/2 − xi−1/2 ,
(3.23)
and the flux limiting function Φ according to van Leer (1985) is defined as Φ(ri+ ) =
|ri+ | + ri+ . 1 + |ri+ |
(3.24)
The argument ri+ of the function Φ is the so-called upwind ratio of two consecutive solution gradients ri+ =
f˜ i − f˜ i−1 + ε . f˜ i+1 − f˜ i + ε
(3.25)
Analogously, one can formulate the flux (Gf˜ )i−1/2 at the left boundary of the control volume Ωi . The expression has to be evaluated with a small parameter, for example ε = 10−10 , in order to avoid division by zero. There are several other limiting functions proposed in the literature, namely, minmod, superbee and MC limiters, and so forth. Each of these limiters leads to a different high-resolution scheme, see LeVeque (2002). The above flux-limited formulation for the scheme guarantees the monotonicity and avoids the possible occurrence of wiggles and negative solution values. In case of discontinuous solutions the scheme gives second-order accuracy in smooth regions while first-order accuracy in the vicinity of shock discontinuities. This formula cannot be applied up to and including boundaries. Therefore, one has to use the first-order upwind formulation in the first two left boundary cells and in the last cell at the right boundary. Again, the resultant ODEs can be solved with any standard ODE-solver, for example the RK45 method in our case which is an embedded Runge-Kutta method of order four and five.
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
4. Numerical test problems The numerical schemes derived in this article can be applied to solve various combinations of nucleation, growth and aggregation processes for different choices of nucleation, growth and aggregation rates. However, to validate the current schemes, it is needed to test the schemes for the test problems where analytical solutions are also available. In the literature, analytical solutions of the PBEs are only available for the following combinations of processes: • Pure growth. • Simultaneous growth and aggregation. • Simultaneous growth and nucleation. In this section, we will consider test cases for above combinations of processes. Note that, analytical solutions for pure aggregation also do exist. However, this case was thoroughly discussed in Filbet and Laurenc¸ot (2004) and Qamar and Warnecke (2006) and will be omitted here. In all numerical test problems the lines represent the analytical solutions while symbols are used to represent the numerical solutions. 4.1. Pure growth 4.1.1. Constant growth problem 1 Here we consider constant growth with G(x) = G0 and an exponential initial distribution of the form: N0 x , (4.1) exp − f (0, x) := f0 (x) = x0 x0 where N0 and x0 are the total number of particles and the mean volume of particles at time t = 0, respectively. In both methods we use a geometric grid (3.3) with q = 3. The analytical solution is given by f (t, x) = f0 (x − G0 t).
(4.2)
Here, we take G0 = 1, N0 = 5, x0 = 0.01 and 60 mesh points. In Fig. 1, the numerical results of both methods are compared with the analytical solution. Note that the first plot which shows the MOC results uses a log scale on both axes so that the solutions are not merely shifted in the graph. The MOC results are in excellent agreement with the analytical solution and are hardly distinguishable. The stability and numerical diffusion of the numerical solution which marred the performance of the finite difference techniques are completely absent in the present scheme. The numerical results of FVS are also in agreement with the analytical solution. However, one can see a clear numerical dissipation in the solution. A comparison of the MOC and FVS shows that FVS results are little bit diffusive. Instead of any other scheme which uses some discretization technique for the advection term, the MOC results show a clear advantage where the advection term disappears from the governing equation. The first plot in Fig. 1 shows that the rapid growth rate results in very rapid growth of small particles, while the large particles stay nearly unchanged.
Fig. 1. A comparison of the results for constant growth problem 1.
1581
1582
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
Fig. 2. A comparison of the results for constant growth problem 2.
4.1.2. Constant growth problem 2 We consider a constant growth problem with initial data of the form: ⎧ 1 ⎪ ⎪ √ exp(−500(0.1x−0.3)2 ) if 2< x ≤4, ⎪ ⎪ ⎪ 0.32π ⎪ ⎪ ⎪ ⎨1 if 6< x ≤ 8, f (0, x)= 1−|x−11| if 10< x ≤ 12, ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 1−100(0.1x−1.5)2 if 14< x ≤ 16, ⎪ ⎪ ⎩ 0 elsewhere. We take the growth rate G(x) = 0.1 and divide the volume xmax = 20 into 100 equal subintervals. In Fig. 2, the numerical results of both methods are compared with the analytical solution and a first-order finite volume scheme. The MOC results are in excellent agreement with the analytical solution and are hardly distinguishable. The numerical results of FVS are also in agreement with the analytical solution. However, one can see a clear numerical dissipation in the solution. Finally, the firstorder finite volume scheme results are more diffusive than both schemes. 4.1.3. Linear growth Here G(x) = G0 x with exponential initial distribution (4.1). For the numerical simulation we take G0 = 1, x0 = 0.01 and N0 = 5. The analytical solution is given by Kumar and Ramkrishna (1997): f (t, x) = f0 (x e−G0 t ) e−G0 t .
(4.3)
In both methods we use a geometric grid (3.3) with q = 4 and 200 mesh points. Fig. 3 shows a comparison of numerical and analytical results. Once again the numerical results of MOC are in excellent agreement with the analytical results. However, one can see a visible under estimation in the results of the FVS on right ends of the number density profiles at t = 4 and 10. The numerical results of the finite volume scheme also show numerical dissipation on the left ends of the number density profiles. The numerical results of the first scheme again justifies the use of MOC for growth term.
Fig. 3. A comparison of analytical and numerical results for linear growth problem.
4.2. Simultaneous growth and aggregation In this case, five analytical solutions are available from Ramabhadran, Peterson, and Seinfeld (1976). The cases represent various combinations of constant and linear growth rates, constant and sum aggregation kernels, as well as exponential and gamma initial distributions. Kumar and Ramkrishna (1997) have also considered all these test problems. The specific choices for each case are given in Table 1 and their corresponding analytical solutions are given in Appendix A. Here we have tested our schemes for all five cases. In each case, we take N0 = 5, x0 = 0.01. A geometric grid with q = 6 is chosen to cover the size range of interest at initial time. Simulation results are presented for different simulation times. In both methods we use 100 mesh points for all five test cases. In the following we discuss all these case one by one. 4.2.1. Constant growth and constant aggregation Here we consider G(x) = G0 and β(x, x ) = β0 . The simulation results for cases 1 and 2 of Table 1 are shown in Figs. 4 and 5, respectively. Both figures indicate that the numerical results of
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
1583
Table 1 Various combinations of growth function, aggregation kernel and initial conditions tested for simultaneous growth and aggregation Case
G(x)
β(x, x )
1
1
100
2
1
100
3
x
10
4
x
10
5
x
x + x
f0 (x)
N0 x exp − x0 x0 x N0 x exp − x0 x0 x x 0 N0 exp − x0 x 0 x N0 x exp − x0 x0 x 0 x N0 exp − x0 x0
the first and second methods are in very good agreement with the analytical results and those in Kumar and Ramkrishna (1997) across several order of magnitude. Even at very small values the results looks very promising. However, the numerical results of
Fig. 5. A comparison of analytical and numerical results for simultaneous growth and aggregation, case 2 (Table 1).
Fig. 4. A comparison of analytical and numerical results for simultaneous growth and aggregation, case 1 (Table 1).
the second method are dissipative at the left discontinuity of the number density profile and is in good agreement in the other regions. There are some differences between first scheme results and analytical and numerical results on the far left side of the number density profile for t = 0.01 in Fig. 5 and similar but not visible differences in all other simulations in Fig. 4 are because the analytical solution in this range breaks down, see also Kumar and Ramkrishna (1997). Ramabhadran et al. (1976) have pointed out that for approximately equal growth and aggregation rates, their analytical solution for G(x) = G0 and β(x, x ) = β0 is valid for only large particle sizes. For the aggregation-dominated situation, the range of validity is increased, where as for the growth-dominated situation, the solution breaks down completely. The simulations given in Figs. 4 and 5 are somehow aggregation-dominated cases in order to allow analytical solutions to hold in a reasonably wide range. For this problem
1584
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
Fig. 7. A comparison of analytical and numerical results for simultaneous growth and aggregation, case 4 (Table 1).
Fig. 6. A comparison of analytical and numerical results for simultaneous growth and aggregation, case 3 (Table 1).
with exponential initial distribution the extent of aggregation M0 (t)/M0 (0), as growth processes conserve total number, is 0.288 and 3.998 × 10−4 at times t = 10−2 and 10, respectively. Similarly the extent of growth processes M1 (t)/M1 (0), as the aggregation process conserve mass, is 1.501 and 4.13 at times t = 10−2 and 10, respectively. 4.2.2. Linear growth and constant aggregation In this case we consider G(x) = G0 x and β(x, x ) = β0 . Figs. 6 and 7 show evolving size distributions for cases 3 and 4 of Table 1. The numerical results are once again in very good agreement with the analytical results and those in Kumar and Ramkrishna (1997). Again, there is a visible smearing in the second method results at far left end of the number density profiles, while are in good agreement with analytical and MOC solutions
in other regions. In case of exponential initial distribution the extent of aggregation M0 (t)/M0 (0) is 0.039 and 3.98 × 10−3 at times t = 1 and 10, respectively. Similarly the extent of growth processes M1 (t)/M1 (0) is 2.72 and 2.2 × 104 at times t = 1 and 10, respectively. 4.2.3. Linear growth and sum aggregation Here we consider G(x) = G0 x and β(x, x ) = β0 (x + x ). The numerical results for the case 5 of Table 1 are presented in Fig. 8. It shows that even with size dependent rate functions, the analytical and the numerical results of first and second methods are in excellent agreement. However, in the results of the second method there is a visible smearing at the left end of the number density profiles and very slight over prediction at the right end. The above test simulations clearly show that the proposed techniques are quite robust and yield very good results for simultaneous growth and aggregation. Here the extent of aggregation M0 (t)/M0 (0) is 0.92 and 0.39 at times t = 1 and 3, respectively. Similarly extent of growth processes
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
M1 (t)/M1 (0) is 2.72 and 20.09 at times t = 1 and 3, respectively. 4.3. Simultaneous nucleation and growth In this section we consider some numerical test problems in order to test the proposed schemes for simultaneous nucleation and growth processes. Most of these numerical problems were also considered by Kumar and Ramkrishna (1997) in order to test their scheme. In both methods we use again geometric grid with q = 6. We take 100 mesh points in both numerical examples.
The volume and time ranges are given by 0 ≤ x ≤ 2.0 and 0 ≤ t ≤ 0.5, respectively. The square step initial condition for the number density is given as (Lim et al., 2002): 100 if 0.4 < x ≤ 0.6, f (0, x) = 0.01 elsewhere. Here we consider a constant growth with G = 1.0 and a linear grid of 200 mesh points. The analytical solution is given as f (t, x) =
4.3.1. Constant growth stiff nucleation This problem was considered by Lim et al. (2002). Suppose that the stiff nucleation takes place at a minimum volume as a function of time: f (0, t) = 100 + 106 exp(−104 (t − 0.215)2 ).
1585
⎧ 2 ⎨ 100 − 106 exp(−104 (Gt − x − 0.215) ) if 0.0 < x ≤ Gt, ⎩
100
if 0.4 < x − Gt ≤ 0.6,
0.01
elsewhere.
The results are shown in Fig. 9. The numerical results of both methods are compared with the analytical solutions and firstorder finite volume scheme. The MOC results are in excellent agreement with the analytical solution and are hardly distinguishable. The numerical results of FVS are diffusive as compared to MOC. Finally, the first-order finite volume scheme results are more diffusive than both schemes. 4.3.2. Constant growth and exponential nucleation In this case we take G(x) = G0 , and S(x) = B0 /x0,n exp(−x/x0,n ). As an initial condition we take the exponential initial distribution of the form f0 (x) = N0 /x0 exp(−x/x0 ). For numerical simulation we take G0 = 1, N0 = 10, x0 = 0.01, x0,n = 0.001, and B0 = 105 . The analytical solution for this case is given as (Kumar and Ramkrishna, 1997): xlow x B0 exp − − exp − , f (t, x)=f0 (x − Gt) + G0 x0,n x0,n (4.4)
Fig. 8. A comparison of analytical and numerical results for simultaneous growth and aggregation, case 5 (Table 1).
Fig. 9. A comparison of exact and numerical results for simultaneous nucleation and growth.
1586
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
Fig. 11. A comparison of the results for linear growth and exponential nucleation. Fig. 10. A comparison of the results for constant growth and exponential nucleation.
and Ramkrishna, 1997): f (t, x) = f0 (x e−G0 t ) exp(−G0 t) +
xlow = max(x1 , x − G0 t).
(4.5)
Here x1 represent the center of the first cell in the mesh. Since now nucleation takes place, we have to add a new cell of nuclei size at each time step. In order to keep the number of cells fixed we delete the end cell at the same time level. The numerical results are shown in Fig. 10. Again the numerical results of both methods are in very good agreement with the analytical solutions. The results of the second method are also excellent with a minor dissipation in the numerical result at t = 10−2 .
− exp −
x x0,n
,
xlow B0 exp − G0 x x0,n
xlow = max (x1 , x e−G0 t ) . (4.6)
The numerical results are shown in Fig. 11. The results of the MOC are in excellent agreement with the analytical solution. However, there is visible dissipation in the numerical results of the FVS. From this and previous results one can conclude that MOC become more important when the growth term is not constant. 4.4. Simultaneous nucleation, growth and aggregation
4.3.3. Linear growth and exponential nucleation Here the initial data and the nucleation term have the same formulation as given in the last example but now G(x) = G0 x. Here we take G0 = 1, N0 = 10, x0 = 0.01, x0,n = 0.001, and B0 = 10. The analytical solution for this case is given as (Kumar
Unfortunately, there is no analytical solution available in this case. We consider the constant growth rate G(x) = G0 , exponential nucleation rate S(x) = B0 /x0,n exp(−x/x0,n ) and a constant aggregation kernel β(x, x ) = β0 . As an initial con-
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
1587
5. Conclusions
Fig. 12. A comparison of FVS and MOC + FVS results with each other and analytical results for nucleation, growth and agglomeration.
dition we take the exponential initial distribution of the form f0 (x) = N0 /x0 exp(−x/x0 ). For numerical simulation we take G0 = 1, N0 = 10, x0 = 0.01, B0 = 1, x0,n = 0.001, β0 = 100, and 100 mesh points. The numerical results for both methods are compared with each other in Fig. 12. The second plot in Fig. 12 represents the results of both methods in the absence of aggregation. The first plot shows that first method resolve the discontinuities very well as compared to the second method. A comparison of both plots, that is with and without aggregation, shows that aggregation significantly effects the distribution of particles present initially and those born afterwards. Particularly, the distribution of the nucleated particles, which lie on the left hand side of discontinuity in the size range, is changed completely in the presence of aggregation. In the second plot since the analytical solutions and MOC solutions were completely overlapping, therefore we omit the analytical solution here. The numerical results of both methods are in good agreement with each other. However, the results of the first method are more resolved as compared to the second method where dissipation is visible at the steep discontinuities.
In this article two numerical methods are proposed for the solution of PBEs for simultaneous nucleation, growth and aggregation processes. The first method combines the finite volume scheme (FVS) for aggregation process with a method of characteristics (MOC) for growth process. For handling nucleation term in this method, a cell of nuclei size is added at a given time level. The second method purely uses a finite volume scheme for nucleation, growth and aggregation of particles. Both methods use the same finite volume scheme for the aggregation process. From the numerical test problems, it was found that the use of MOC for linear growth processes is more successful as compared to any finite difference and even finite volume schemes. Although the growth process looks simpler than the aggregation process, still it give rise to complexity in the PBE. Since the presence of growth term introduces hyperbolicity in the governing equation, there may exist discontinuities or steep fronts in the solution. The approximation of advection term for the growth process by a finite difference or a finite volume scheme introduces an inherent viscosity in the solution. This leads to the smearing of the resultant solution in the vicinity of discontinuities (steep fronts). In case of MOC, the solution moves along the propagation path-line and the advection term ∂Gf˜ /∂x simply disappears from the governing equation. Therefore, no special care is needed for the advection term which is usually required in other schemes. Even though the results of MOC are superior to the FVS, still the FVS worked very well for all numerical test cases presented in this article. The use of FVS for the current integrodifferential population balance equation show their applicability, generality and robustness. The schemes also show their flexibility to be adapted in the same manner and to the same grids which were basically proposed for the methods specific to PBEs. Note that, the aggregation term in both schemes is calculated with the same finite volume scheme. This shows that finite volume scheme works very well for aggregation terms and gives very accurate results as for as the advection terms is approximated very well. Since, in the first method the growth term is calculated with the MOC, which gives more accurate results, the overall results for simultaneous growth and aggregation are also very excellent and almost indistinguishable from the analytical results. Especially, for the two-component problems involving simultaneous growth and aggregation processes, the use of pure finite volume seems to be computationally expensive. It was found that, the use of MOC for growth terms makes the first method at least two times faster than the second method. However, in most of the test cases presented in this article the computational time was not more than a minute for both methods. There were only one or two examples where the computational time was a little bit more than a minute. Several numerical results in different combinations of nucleation, growth and aggregation are presented in this article. The numerical results of the schemes were found to be in good agreement with the analytical solutions.
1588
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589
Here Λ = G0 /β0 N0 x0 and Q + Q2 − R g1 = − , Rx0
Table 2 Moments M0 and M1 of the distributions Case
Moments
β(x, x ) = β0
M0 =
G(x) = G0
M1 = N0 x0 1 −
β(x, x ) = β0
M0 =
G(x) = G0 x β(x, x )
= β0 (x
2N0 2 + β0 N0 t
2G0 ln β0 N0 x0
2 2 + β0 N0 t
2N0 2 + β0 N0 t
M1 = N0 x0 exp(β0 t) + x )
G(x) = G0 x
M0 = N0 exp
β N x 0 0 0
G0 M1 = N0 x0 exp(β0 t)
(1 − exp(G0 t))
g2 = −
Q(x) = χ − Λ(χ − 1 − χ ln χ),
Q−
χ=
Q2 − R , Rx0
N0 , M0
R(χ) = χ − (2Λ)2 χ(χ − 1 − ln χ) − (χ − 1)2 .
(A.5)
A.2. Linear growth and constant aggregation In case of linear growth G(x) = G0 x and size-independent aggregation kernel β(x, x ) = β0 , the analytical solution for exponential initial distribution (A.1) is given by M02 M0 f (t, x) = exp − x . (A.6) M1 M1
Appendix A. Analytical solutions for simultaneous growth and aggregation problems There are several combinations of growth rates and aggregation kernels where analytical solutions are available. These solution were derived by Ramabhadran et al. (1976). Since two simultaneous processes are taking place, the dynamics of the particle size distribution depends on two characteristic times, one for aggregation and other for growth. The ratio of these two times is the basic dimensionless parameter of the problem. Let us define Λ as the ratio of the characteristic times of growth and aggregation. The analytical solution for two initial distribution and different combinations of the growth and aggregation rates are given below.
For Gaussian-type initial distribution (A.3) the analytical solution is given as f (t, x) =
M02 1 √ M1 1 − M0 /N0 N0 x M 0 N0 x . sinh × exp − 1− M1 N0 M1
(A.7)
The moments M0 and M1 for this case are given in Table 2. A.3. Linear growth and sum aggregation
A.1. Constant growth and constant aggregation For the constant growth G(x) = G0 and size-independent aggregation kernel β(x, x ) = β0 , the analytical solution with exponential initial distribution: N0 x , (A.1) exp − f (0, x) := f0 (x) = x0 x0 is given as f (t, x) =
M02 /M1 1 − 2Λx0 (N0 − M0 /M1 ) M0 (x − 2Λx0 (N0 /M0 − 1)) × exp − , M1 1 − 2Λx0 (N0 − M0 /M1 )
For the case that G(x) = G0 x and β(x, x ) = β0 (x + x ), the analytical solution for exponential initial distribution (A.1) is given as f (t, x) =
M0 √ 1 − M0 /N0 M0 2N0 M0 N0 x , × exp − −1 x I1 2 1− M1 M0 N0 M1 x
(A.8) (A.2)
where Λ = G0 /β0 N0 x0 . The moments M0 and M1 are given in Table 2. For Gaussian-type initial distribution N0 x x , (A.3) f (0, x) := f0 (x) = 2 exp − x0 x0 the analytical solution is given by M0 /x02 N0 f (t, x) = −1 exp g1 x − 2Λx0 R(g1 − g2 ) M0 N0 − exp g2 x − 2Λx0 −1 . (A.4) M0
where I1 is the modified Bessel function of first kind of order one. The moments M0 and M1 for are given in Table 2. References Abramowitz, M., & Stegun, I.A. (1972). Handbook of mathematical functions with formulas, graphs and mathematical tables (p. 11). New York: Dover. David, R., Villermaux, P., & Klein, J. P. (1991). Crystallization and precipitation engineering-IV. Kinetic model to adipic acid crystallization. Chemical Engineering Science, 46, 1129–1136. Filbet, F., & Laurenc¸ot, P. (2004). Numerical simulation of the Smoluchowski coagulation equation. SIAM Journal on Scientific Computing, 25(6), 2004–2048. Gunawan, R., Fusman, I., & Braatz, R. D. (2004). High resolution algorithms for multidimensional population balance equations. AIChE Journal, 50, 2738–2749. Hounslow, M. J. (1990). A discretized population balance for continuous systems at steady-state. AIChE Journal, 36, 106–116.
S. Qamar, G. Warnecke / Computers and Chemical Engineering 31 (2007) 1576–1589 Hounslow, M. J., Ryall, R. L., & Marshall, V. R. (1988). A discretized population balance for nucleation, growth, and aggregation. AIChE Journal, 34, 1821–1832. 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. Lapidus, L., & Pinder, G. F. (1982). Numerical solutions of partial differential equations in science and engineering. New York, USA: Wiley. LeVeque, R. J. (2002). Finite volume methods for hyperbolic problems. New York, NY: Cambridge University Press. 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. Litster, J. D., Smith, D. J., & Hounslow, M. J. (1995). Adjustable discretized population balance for growth and agglomeration. AIChE Journal, 41, 591–603. 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 and Chemical Engineering, 26, 1103–1116. Motz, S., Mitrovi´c, A., & Gilles, E. D. (2002). Comparison of numerical methods for the simulation of dispersed phase systems. Chemical Engineering Science, 57, 4329–4344. Muhr, H., David, J., Villermaux, J., & Jezequel, P. H. (1996). Crystallization and precipitation engineering VI: Solving population balance in the case of the precipitation of silver bromide crystals with high primary nucleation rate by
1589
using first order upwind differentiation. Chemical Engineering Science, 51, 309–319. Qamar, S., Ashfaq, A., Warnecke, G., Angelov, I., Elsner, M. P., & SeidelMorgenstern, A. (2007). Adaptive high resolution schemes for population balances in crystallization processes. Computers and Chemical Engineering, 31, 1296–1311. 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 and Chemical Engineering, 30, 1119–1131. Qamar, S., & Warnecke, G. (2006). Solving population balance equations for two-component aggregation by a finite volume scheme. Chemical Engineering Science, 62, 679–693. Ramabhadran, T. E., Peterson, T. W., & Seinfeld, J. H. (1976). Dynamics of aerosol coagulation and condensation. AIChE Journal, 22, 840– 851. Ramkrishna, D. (1985). The status of population balances. Reviews in Chemical Engineering, 3(1), 49–95. 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. 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. van Leer, B. (1985). Upwind-difference methods for aerodynamic problems governed by the Euler equations. Lectures in applied mathematics (Vol. 22-Part 2, pp. 327–336). Providence, RI: American Mathematical Society.
