On the solution of population balances for nucleation, growth, aggregation and breakage processes

Chemical Engineering Science 64 (2009) 2088 -- 2095

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / c e s

On the solution of population balances for nucleation, growth, aggregation and breakage processes Shamsul Qamar a,∗ , Gerald Warnecke b , Martin Peter Elsner c a b c

Department of Mathematics, COMSATS Institute of Information Technology, Plot 30, H-8/1 Islamabad, Pakistan Institute for Analysis and Numerics, Otto-von-Guericke University Magdeburg, 39106 Magdeburg, Germany Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstrasse 1, 39106 Magdeburg, Germany

A R T I C L E

I N F O

Article history: Received 15 August 2008 Received in revised form 1 December 2008 Accepted 19 January 2009 Available online 4 February 2009 Keywords: Population balance models High resolution finite volume schemes Method of characteristics Nucleation Growth Aggregation Breakage Hyperbolic conservation laws

A B S T R A C T

This work is concerned with the modeling and simulation of population balance equations (PBEs) for combined particulate processes. In this study a PBE with simultaneous nucleation, growth, aggregation and breakage processes is considered. In order to apply the finite volume schemes (FVS) a reformulation of the original PBE is introduced. This reformulation not only help us to treat the aggregation and breakage processes in a manner similar to the growth process in the FVS but also in deriving a stable numerical scheme. Two numerical methods are proposed for the numerical approximation of the resulting reformulated PBE. The first method combines a method of characteristics (MOC) for growth process with an FVS for aggregation and breakage processes. The second method purely uses a semidiscrete FVS for all processes. Both schemes use the same FVS for aggregation and breakage processes. The numerical results of the schemes are compared with each other and with the available analytical solutions. The numerical results were found to be in good agreement with analytical solutions. © 2009 Elsevier Ltd. All rights reserved.

1. Introduction Population balance models (PBMs) are encountered in several scientific and engineering disciplines describing the temporal change of the number density function. These integro-differential equations describe the change due to continuous transport of particles in the state space, e.g. by particle growth as well as the effect of birth and death events such as nucleation, breakage and aggregation processes. They were introduced in the field of chemical engineering by Hulbert and Katz (1964) and later fully developed by Randolph and Larson (1988). However, their applications were limited due to a lack of computational power. Since high speed computers are available now, they became popular which is reflected by their use in several areas of chemical and biochemical engineering. They are used to study precipitation, polymerization, crystallization, food processes, pharmaceutical manufacture, pollutant formation in flames, particle size distribution (PSD) of crushed material and rain drops, dispersed phase distributions in multiphase flows, and growth of microbial and cell populations. In practical situations the population density function may extend over orders of magnitudes and can be very sharp. Therefore,

∗ Corresponding author. Tel.: +92 51 9258481x814; fax: +92 51 4442805. E-mail address: [email protected] (S. Qamar). 0009-2509/$ - see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2009.01.040

accurate numerical simulation of the population density functions can be challenging for a given numerical scheme. Due to this fact several researchers working in this field have developed specialized algorithms for solving population balance equations (PBEs), see for example Ramkrishna (1985, 2000) and Nicmanis and Hounslow (1998), and references therein. As a result, several numerical methods were developed such as Monte Carlo methods by Lee and Matsoukas (2000) and Mishra (2000), the finite differences/ discretized population balances by Hill and Ng (1995), Hounslow et al. (2001), Kumar et al. (2006), Kumar and Ramkrishna (1996) and Vanni (1999), the quadrature method of moments by Daniele et al. (2003), Kostoglou and Karabelas (2002, 2004) and the finite element methods by Everson et al. (1997). Kumar and Ramkrishna (1997) have combined their discretization technique on a non-uniform grid for aggregation terms with the method of characteristics (MOC) for growth terms in order to solve the above PBEs for simultaneous nucleation, growth and aggregation processes. While, Kumar et al. (2006) have used the cell averaged technique, which is an improved version of the fixed pivot technique, to solve different combinations of the above processes. Moreover, Lim et al. (2002) have combined the same discretization technique of Kumar and Ramkrishna (1997) with the weighted essentially non-oscillatory (WENO) scheme as well as with the MOC for the simulation of potassium sulfate problem (K2 SO4 –H2 O system) involving nucleation, growth, aggregation and breakage

S. Qamar et al. / Chemical Engineering Science 64 (2009) 2088 -- 2095

kinetics. They have found that MOC gives better results than the WENO scheme for the approximation of growth process. Filbet and Laurençot (2004) have proposed a conservative finite volume schemes (FVS) for single-component aggregation problem. These authors have rewritten the PBE for aggregation problems in a form which can be readily solved by an FVS. Later on, Qamar and Warnecke (2006) have extended the scheme in order to simulate two-component aggregation problems. High resolution FVS, originally developed for gas dynamics, were also used for the simulation of PBMs describing crystallization process, see Ma et al. (2002), Gunawan et al. (2004), and Qamar et al. (2006, 2007). Recently, Gunawan et al. (2008) have presented the parallel implementation of the high resolution FVS for simulating general PBMs incorporating nucleation, growth, aggregation and breakage processes. Moreover, Woo et al. (2006) have simulated the mixing effect in antisolvent crystallization by coupling turbulent computational fluid dynamics (CFD) code with a multi-environment probability density function (PDF) model for capturing the micromixing in the subgrid scale, and with the PBE for the evolution of crystal size distribution. They have used the semidiscrete central for the discretization of PBE along the internal coordinate. The proposed high resolution finite scheme presented in this article can also be used for such time of coupling of PBE with CFD code in order to account for external variations. We expect better results from the current high resolution schemes approximating the growth process due to its third order accuracy, see Koren (1993). The aim of this current work is to numerically solve PBEs containing nucleation, growth, aggregation and breakage processes. Two methods are used to solve the given PBE. The first method combines a MOC for growth term, see Kumar and Ramkrishna (1997), with an FVS for aggregation and breakage terms. In the second method, a semi-discrete FVS is used for all processes, see also Note that, in both methods the aggregation kinetics is solved with the same finite volume formulation, see Qamar and Warnecke (2007). Themain difference in both methods is the way they calculate the growth term. The first method uses MOC for the growth term, while the second method uses an FVS of Koren (1993) to discretize the growth term. For handling the nucleation term in the first method, MOC is combined with a procedure of adding a cell of the nuclei size at each time level. A standard ordinary differential equation (ODE) solver can be used to solve the resultant ODEs. 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 processes. Section 4, presents numerical results. Finally, Section 5 gives conclusions and remarks.

2. Mathematical model The PBE for an ideally mixed system describes the dynamics of the volume distribution function f := f (t, x)  0 of particles of volume x > 0 at time t  0. The one-dimensional PBE for a well mixed system is given as Perry and Green (1997)

jf (t, x) j[(Gf )(t, x)] ± + = Q± agg (t, x) + Qbreak (t, x) + Qnuc (t, x), jt jx f (0, x) = f0 (x),

x ∈ R+ ,

2089

± are nucleated. The aggregation term Q± agg := Qagg (t, x) is defined as

Q± agg =



1 2 −

x

(t, x − x , x )f (t, x − x )f (t, x ) dx

0 ∞ 0

(t, x, x )f (t, x)f (t, x ) dx .

(3)

The first term in the above equation 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 volume 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 0  (t, x, x ) = (t, x , x),

(x, x ) ∈ R2+ .

It is basically the probability that particles of sizes x, x meet and stick together. ± Lastly the breakage term Q± break := Qbreak (t, x) is given as  ∞  Q± b(t, x, x )S(x )f (t, x ) dx − S(x)f (t, x). (4) break = x

The breakage function b(t, x, x ) is the PDF for the formation of particles of size x from particle of size x . The selection function S(x ) describes the rate at which particles are selected to break. The breakage function has the following properties:  x  x   ˜ b(t, x, x ) dx = N(x) and y b(t, y , x) dy = x, (5) 0

0

˜ where the function N(x) represents the number of fragments obtained from the breakage of a particle of volume x. The moments of f are defined as  ∞ xi f (t, x) dx, i = 0, 1, 2, . . . , (6) Mi (t) = 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. To apply an FVS, it is convenient to reformulate PBE (1) in a mass balance form. If we apply the FVS directly to (1) then we have to consider all the terms on the right-hand side as cell averaged source terms. This may result in an unstable numerical scheme if the aggregation and breakage processes are dominant compared to the growth process. To avoid this problem, the idea is to treat the aggregation and breakage processes as cell interface fluxes similar to the growth process in the finite volume. Therefore, the resulting reformulation not only help us to treat the aggregation and breakage processes in a manner similar to the growth process in the FVS but also in deriving a stable numerical scheme. Multiplying both sides of Eq. (1) by x, we get x

jf (t, x) j[(Gf )(t, x)] ± +x = xQ± agg + xQbreak + xQnuc , jt jx x ∈ R+ .

xf (0, x) = xf 0 (x),

(t, x) ∈ R2+ , (7) (8)

(1)

˜ Let us define f˜ (t, x) := xf (t, x) and Q nuc (t, x) := xQnuc (t, x). The aggregation term xQagg (t, x) can be rewritten as Filbet and Laurençot (2004), Qamar and Warnecke (2007)

(2)

xQagg (t, x) = −

where (t, x) ∈ R2+ and R+ :=]0, +∞[. Here f (t, x) represents the number density per unit volume of the homogeneous dispersion, G(t, x) is the growth term, Qnuc (t, x) is the rate at which particle of size x

where

Fagg (t, x) = −

jFagg (t, x) , jx 

x 0



∞

x−u

u(t, u, v)f (t, u)f (t, v) dv du.

(9)

(10)

2090

S. Qamar et al. / Chemical Engineering Science 64 (2009) 2088 -- 2095

Similarly, the breakage term according to Appendix A has the form xQbreak (t, x) = where

Fbreak (t, x) =

jFbreak (t, x) , jx 

x

0



∞

x

(11)





y b(t, y , x )S(x )f (t, x ) dx dy .

(12)

The same Eqs. (11) and (12) were also obtained by Kumar (2006) with a different approach. If x is independent of time t then Eq. (7) in the light of above definitions gives

jf˜ (t, x) j[(Gf˜ )(t, x)] G(t, x)f˜ (t, x) + − jt jx x jFagg (t, x) jFbreak (t, x) ˜ + + Qnuc (t, x) =− jx jx

(13)

with initial data f˜ (0, x) = f˜0 .

(14)

However, if x is a function of time, i.e. x := x(t), we obtain x

jf (t, x) jf˜ (t, x) dx = + f (t, x) jt jt dt =

jf˜ (t, x) G(t, x)f˜ (t, x) + . jt x

(15)

Hence, the corresponding reformulation of Eq. (7) gives

jFagg (t, x) jFbreak (t, x) ˜ jf˜ (t, x) j[(Gf˜ )(t, x)] + =− + +Qnuc (t, x). jt jx jx jx

(16)

3. Numerical procedure In the following we will derive two methods for solving the above modified PBEs (13) and (16). 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, i.e. regular (or irregular grid) and geometric grid. Regular/irregular grid: Let N be a large integer representing the total number of cells, we denote by (xi−1/2 ) the partitions i∈{1, ...,N+1} of [xmin , xmax ], where xmin is the minimum and xmax is the maximum value of the property variable. We set for i = 1, 2, . . . N − 1 xN+1/2 = xmax ,

xi+1/2 = xmin + ixi .

(17)

Furthermore, we set xi = (xi−1/2 + xi+1/2 )/2,

xi = xi+1/2 − xi−1/2 .

(18)

Geometric grid: In order to capture the initial profile properly, in some test problems it is useful to use a geometric grid. In other words, geometric grid is useful when smaller particle volume range is more important compared to the larger particle range. A typical one-dimensional geometric grid is given as x1/2 = xmin ,

xi+1/2 = xmin + 2(i−N)/q (xmax − xmin ),

 1 f˜i (0) = f˜ (x) dx. xi i 0

(20)

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 (13) and (14). 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 MOC for the growth process and an FVS for the aggregation and breakage processes. 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 pathline, the advection term jGf˜ / jx 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 (17) or (19) 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

Instead of the original Eq. (1) we will discretize (13) or (16) in our numerical schemes. Due to the relation f˜ (t, x) := xf (t, x) one can easily recover f (t, x).

x1/2 = xmin ,

Let i = [xi−1/2 , xi+1/2 ] for i  0. We approximate the initial data f˜0 (x) in each grid cell by

i = 1, 2, . . . , N. (19)

The value of parameter q describes the extent of grid refinement in a smaller particle range.

(21)

˜ Hence the reformulation of Eq. (16) for f˜ := f˜ (t, x) and Q nuc := ˜ Qnuc (t, x) is given as 



jFagg jFbreak ˜ jf˜ j dx ˜ + + + Qnuc . f =− jt jx dt jx jx

(22)

˜ ˜ Let f˜i = f˜i (t) and (Q nuc )i = (Qnuc )i (t) denote, respectively, the average values of the number density and nucleation term in each cell i as f˜i :=

 1 f˜ dx, xi (t) i (t)

˜ (Q nuc )i :=

 1 Q˜ dx. xi (t) i (t)

(23)

Using the Leibniz integration rule (Abramowitz and Stegun, 1972) and the above definitions, integration of Eq. (22) over the control volume i (t) = [xi−1/2 (t), xi+1/2 (t)] gives after simplification (Qamar and Warnecke, 2007) 1 df˜i = − [(Fagg )i+(1/2) − (Fagg )i−(1/2) ] dt xi (t) 1 [(Fbreak )i+(1/2) − (Fbreak )i−(1/2) ] + xi (t) f˜ ˜ − (Gi+(1/2) − Gi−(1/2) ) i + (Q nuc )i , xi (t)

(24)

where the cell interface values (Fagg )i+(1/2) according to Filbet and Laurençot (2004) are given as (Fagg )i+1/2 =

i  k=0

xk (t)f˜k

⎧ N ⎨ ⎩

j=i,k

f˜j

 hj

(x , xk ) x

dx



⎫  i,k −1/2 ⎬ , x )  (x  k + f˜i,k −1 dx + O(x3i ).  ⎭ x xi+1/2 −xk

(25)

Here, the integer i,k corresponds to the index of the cell such that xi+1/2 (t) − xk (t) ∈ i,k −1 (t).

S. Qamar et al. / Chemical Engineering Science 64 (2009) 2088 -- 2095

Similarly the cell interface values (Fbreak )i+(1/2) are given by (Fbreak )i+(1/2) =

⎛

i   k=0

k

x∗ ⎝

N 

f˜j



j=i+1

j

⎞

b(x∗ , x )

S(x )  ⎠ ∗ dx dx + O(x3i ). x (26)

In summary, we have to solve the following set of equations: df˜i 1 = − [(Fagg )i+(1/2) − (Fagg )i−(1/2) ] dt xi (t) 1 [(Fbreak )i+(1/2) − (Fbreak )i−(1/2) ] + xi (t) f˜ ˜+ ), − (Gi+(1/2) − Gi−(1/2) ) i + (Q nuc i xi (t) dx˜ i+(1/2) dt

(27)

2091

where xi+(1/2) = xi+1 − xi . The argument ri+(1/2) of the function  is the so-called upwind ratio of two consecutive solution gradients ri+(1/2) =

Fi+1 − Fi +  . Fi − Fi−1 + 

(34)

Analogously, one can formulate the flux Fi−(1/2) at the left boundary of the control volume i . The expression has to be evaluated with a small parameter, e.g.  = 10−10 in our numerical experiments, in order to avoid division by zero. There are several other limiting functions proposed in the literature, namely, minmod, superbee, MC limiters, etc. Each of these limiters leads to a different high resolution scheme, see LeVeque (2002). Again, a standard ODE-solver can be used to solve the resulting system of ODEs. 4. Numerical test problems

= Gi+(1/2) ,

∀i = 1, 2, . . . , N

(28)

with initial condition f˜ (0, xi ) = f˜0 (xi ).

(29)

The above system of ordinary differential equations can be solved by any standard ODE solver. 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 occurs at the minimum particle size, one can use the nucleation simply as a boundary condition xmin B0 (t) . f˜ (t, xmin ) = G(t, xmin )

(30)

Here, B0 (t) is the nucleation rate of particles of minimum particle size xmin . 3.2.2. Method II: semidiscrete FVS Alternatively, one can use a semidiscrete FVS in order to solve PBE (13) without any further modification. ˜ ˜ Let f˜i =f˜i (t) and (Q nuc )i =(Qnuc )i (t) denote, respectively, the average value of the number density and the nucleation term as given in (23). Moreover, let F(t, x) := (Gf )(t, x). Integration of (13) over the control volume i =[xi−(1/2) , xi+(1/2) ] gives Qamar and Warnecke (2006, 2007) 1 df˜i 1 =− [F −Fi−(1/2) ]− [(Fagg )i+(1/2) −(Fagg )i−(1/2) ] dt xi i+(1/2) xi Gi+(1/2) f˜i 1 ˜ [(Fbreak )i+(1/2) −(Fbreak )i−(1/2) ]+ , (31) +(Q nuc )i + xi xi where (Fagg )i+(1/2) and (Fbreak )i+(1/2) are given by (25) and (26), respectively. The flux Fi+(1/2) at the right cell interface is given by the following flux-limited formula (Koren, 1993)

Fi+(1/2) = Fi + (1/2)(ri+ )(Fi − Fi−1 )

(32)

and the flux limiting function  according to Koren (1993) is defined as    1 xi (ri+(1/2) ) = max 0, min 2ri+(1/2) , min 3 xi−(1/2)  2ri+(1/2) xi + ,2 , (33) 3 xi+(1/2)

The numerical schemes derived in this article can be applied to solve various combinations of nucleation, growth, aggregation and breakage processes. In this section, we consider three test problems. 4.1. All processes simultaneously This problem was considered by Lim et al. (2002). Suppose that the stiff nucleation takes place at a minimum particle size (xmin = 0) as a function of time f (t, 0) = 100 + 106 exp(−104 (t − 0.215)2 ).

(35)

Hence, we consider the nucleation as a left boundary condition. The particle size and time ranges are 0  x  2.0 and 0  t  0.5, respectively. The square step initial data for the number density is given as  100 for 0.4  x  0.6, f (0, x) = (36) 0.01 elsewhere. We consider constant growth rate with G = 1.0. The analytical solution for only growth and nucleation is given as Lim et al. (2002) ⎧ 2 6 4 ⎨ 100+10 exp(−10 ((Gt−x)−0.215) ) for 0.0  x  Gt, f (t, x)= 100 for 0.4  x−Gt  0.6, ⎩ 0.01 elsewhere. In this solution, a square step discontinuous shock and a narrow wave which is originated from nucleation move along the propagation path-line, x = x0 + Gt. The numerical test is carried out on 200 grid points. The kinetic parameters from aggregation and breakage are  = 1.5 × 10−5 , b(t, x, x ) = 2/x and S(x) = x2 . Fig. 1, which is obtained from Method I (MOC+FVM), depicts the effects of the growth, nucleation, aggregation, and breakage terms on the PSD. The solid line is the analytic solution for pure growth and nucleation problem but without aggregation and breakage processes. The numerical results of MOC for pure growth and nucleation are overlapping with the analytical solution as shown in Fig. 1. The stiff nucleation at the left boundary, produces a sharp peak. The aggregation term causes peaks of the PSD to be smeared and increases the population of large size particles. In contrast, the breakage term increases the population of small size particles. Therefore, the PSD of the PBE with the four kinetics is dispersed more broadly. Similar results which are obtained from the Method II (FVS) are given in Fig. 2. In comparison to Method I, more numerical diffusion is visible in the steep gradients of the Method II solutions. However, the high resolution scheme still resolves all the profiles of the solution quite well. Table 1 reports the computational time of Methods I and II for different combinations of processes. The table presents that MOC

2092

S. Qamar et al. / Chemical Engineering Science 64 (2009) 2088 -- 2095

Pure Growth and Nucleation

Nucleation + Growth + Aggregation 106

MOC Analytical Pure Growth and Nucleation

104

number density

number density

106

102 100 10−2

MOC+FVM Analytical Pure Growth and Nucleation

104 102 100 10−2

0

0.5

1 volume

1.5

2

0

number density

number density

106

MOC Analytical Pure Growth and Nucleation

104

1 volume

1.5

2

Nucleation + Growth + Aggregation + Breakage

Nucleation+Growth + Breakage 106

0.5

102 100 10−2

MOC+FVM Analytical Pure Growth and Nucleation

104 102 100 10−2

0

0.5

1 volume

1.5

2

0

0.5

1 volume

1.5

2

Fig. 1. Results of Method I (MOC+FVS) at t = 0.5 and N = 200 mesh points.

Pure Growth and Nucleation

Nucleation + Growth + Aggregation 106

FVM Analytical Pure Growth and Nucleation

104

number density

number density

106

102 100 10−2

FVM Analytical Pure Growth and Nucleation

104 102 100 10−2

0

0.5

1 volume

1.5

2

0

number density

number density

106

FVM Analytical Pure Growth and Nucleation

104

1 volume

1.5

2

Nucleation +Growth + Aggregation + Breakage

Nucleation + Growth + Breakage 106

0.5

102 100 10−2

FVM Analytical Pure Growth and Nucleation

104 102 100 10−2

0

0.5

1 volume

1.5

2

0

0.5

1 volume

Fig. 2. Results of Method II (FVS) at t = 0.5 and N = 200 mesh points.

1.5

2

S. Qamar et al. / Chemical Engineering Science 64 (2009) 2088 -- 2095

is faster for the growth processes as compared to FVM. However, the overall computational cost is comparable for both methods after including other processes. The total L1 -error in the case of growth and nucleation processes is 104 for MOC and 1.92 × 104 for FVS. Similarly, the total L1 -error in the case of all processes is 6.63 × 103 for Method I and 4.76×104 for Method II with 100 grid points. Since, there is no analytical solution for the last case, we have taken the solution of Method I for 200 grid points as a reference solution. These numerical errors show that for the same computational time and grid points Method I gives more accurate results compared to Method II. Furthermore, our numerical methods are comparable in accuracy to those presented in Lim et al. (2002) but are faster according to their CPU time calculations. Table 1 Calculation time for Methods I&II. Description

CPU: Method I (s) −10

Growth and nucleation Growth, nucleation and breakage Growth, nucleation and aggregation All processes

10 3.76 20.73 24.56

CPU: Method II (s) 0.1 3.62 20.81 24.51

Table 2 Initial data and analytical solutions for pure breakage process.

2093

4.2. Pure breakage Here we consider two test cases with different initial data that can be solved analytically. The analytical solutions for different initial PSDs and other parameters have been provided by Ziff and McGray (1985). In these test problems we take the maximum particle size Rx = 125. The analytical solutions and the corresponding initial conditions are given in Table 2. In the first test case, exponential distribution is taken as initial PSD. In Fig. 3, the analytical and numerical PSD along with their first three moments are compared with each other. The final simulation time is t = 10. Here we have used the geometric grid with q = 3 and 60 mesh points. The results for the PSD are given on semi-log scales while linear scale is used for the moments plots. Since particles flows towards the small size range in the breakage process, it is very important to see the performance of the scheme there. From the results it is clear that FVS predict the PSD and the first three moments very well. In the second test case, again the exponential distribution is taken as initial PSD. Fig. 4 shows the results. The final simulation time is t = 5. The numerical results are in very good agreement with analytical solutions. One cannot see difference between the analytical and numerical solutions. Again we have used geometric grid with q = 3 and 60 mesh points. 4.3. Simultaneous aggregation and breakage

Case

S(x)

b(x, x )

f (0, x)

Analytical solution, f (t, x)

1 2

x x2

2/x 2/x

exp(−x) exp(−x)

exp(−x(1 + t))(1 + t) exp(−tx2 − x)[1 + 2t(1 + x)]

2

Here, we consider simultaneous aggregation and breakage processes. The numerical results are compared with the available analytical solutions of Lage (2002). We consider a uniform binary

12

140

FVS Analytical

120

FVS Analytical

10 Mj (t)/Mj (0)

number density

j=0 100 80 60 40

8 6 4 2

20

j=1

0

0 10−6

10−4

10−2 volume

100

102

j=2 0

2

4

6

8

10

time

Fig. 3. Pure breakage problem, Case 1 in Table 2, q = 3, N = 60 and t = 10.

4.5

12

FVS Analytical

3.5 Mj (t)/Mj (0)

number density

10

FVS analytical

4

8 6 4

j=0

3 2.5 2 1.5

j=1

1 2 0 10−4

0.5 10−3

10−2

10−1 100 volume

101

102

j=2

0 0

1

2

3 time

Fig. 4. Pure breakage problem, Case 2 in Table 2, q = 3, N = 60 and t = 5.

4

5

2094

S. Qamar et al. / Chemical Engineering Science 64 (2009) 2088 -- 2095

FVS Analytical

100

number density

number density

10−2 10−4 10−6 10−8 10−10 10−3

10−2

10−1 100 volume

101

102

1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

FVS Analytical

10−3

10−2

10-1 100 volume

101

102

Fig. 5. Simultaneous aggregation and breakage with initial condition (37) at t = 6. Left-hand side plot is on log–log scale and the right-hand side on semi-log scale.

Simultaneous Growth, Nucleation and Agglomeration

104

t=0 t = 1e−4 t = 1e−4 t = 1e−3 t = 1e−3 t = 1e−2 t = 1e−2

102 100

Simultaneous Growth and Nucleation t=0 t = 1e−4 t = 1e−4 t = 1e−3 t = 1e−3 t = 1e−2 t = 1e−2

103 102 number density

number density

104

10−2

101 100 10−1 10−2 10−3

10−4

10−4 10−4

10−2 volume

100

102

10−4

10−3 10−2 volume

10−1

100

Fig. 6. A comparison of the results of FVS (symbols) and MOC+FVS (lines) with each other for nucleation, growth and aggregation processes.

breakage b(x, x ) = 2/x , linear selection function S(x) = S0 x, and constant aggregation kernel (x, x ) = 0 . The initial condition is given according to the following distribution function: f (0, x) =

  M02 M0 exp − x . M1 M1

(37)

A geometric grid discretization has been used for this problem with q=3 and N=60. The final simulation time is t=6. Fig. 5 shows the numerical results on both log–log and semi-log scales. In this figure the left-hand side results are given on log–log scale, while on the right-hand side the same results are plotted on semi-log scale. On the log–log scale the FVS seems in agreement with analytical solution. However, the numerical results on semi-log scale show visible under estimation in smaller volume range. In contrast, it was observed by Kumar (2006) that for the same problem the fixed pivot technique gives an over predicting solution in the smaller volume range. On the other hand, the cell averaged technique of Kumar (2006), which is an improved version of the fixed pivot technique, has given almost overlapping results with the exact solutions. 4.4. Simultaneous nucleation, growth and aggregation Unfortunately, there is no analytical solution available in this case. We consider the constant growth rate G(x) = G0 , exponential nucleation rate Q+ nuc (x) = B0 /x0,n exp(−x/x0,n ) and a constant aggregation kernel (x, x ) = 0 . 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, 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. 6. The second plot in Fig. 6 represents the results of both methods in the absence of aggregation process. The first plot shows that Method I (FVS+MOC) resolve the discontinuities very well as compared to Method II (FVS). A comparison of both plots, i.e. with and without aggregation process, 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 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. 5. Conclusions In this article two numerical methods are proposed for the solution of PBEs for simultaneous nucleation, growth, aggregation and breakage processes. The first method combines the FVS for aggregation and breakage processes with the MOC for growth process. The second method purely uses an FVS for all processes. 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 an FVS introduces an inherent viscosity in the solution.

S. Qamar et al. / Chemical Engineering Science 64 (2009) 2088 -- 2095

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 jGf˜ / jx 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 quite well for some of the numerical test cases presented in this article. Note that, aggregation and breakage terms in both schemes are calculated with the same FVS. This shows that FVS works quite well for aggregation and breakage terms and gives accurate results as far as the advection terms is well approximated. 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 excellent. Moreover, the role of FVS will be more important when one combines the PBMs with CFD models. In that case one can use MOC for the growth process and FVS for the spatial advection terms. It was found that for the same number of grid points the Method I gives better solution compared to Method II and the simulation times were comparable. In most of the test cases presented in this article the computational time was not more than a minute for both methods. Finally, it is important to mention that current scheme is also useful for coupled CFD–PBE models. The proposed numerical scheme can be easily coupled with any CFD code in order to simulate nonideally mixed systems. This is also a topic of our future study. Appendix A. Verification of Eq. (11) In the following we show that how one can obtain the left-hand side of Eq. (11) from the right-hand side. Using Eq. (12) in (11) we get ˜ jF j break (t, x) = jx jx 

=

x



∞



x

0

 x

∞



y b(t, y , x )S(x )f (t, x ) dx dy

xb(t, x, x )S(x )f (t, x ) dx x



 

− y b(t, y , x)S(x)f (t, x) dy 0  ∞  = xb(t, x, x )S(x )f (t, x ) dx x  x  − S(x)f (t, x) y b(t, y , x) dy 0  ∞  = xb(t, x, x )S(x )f (t, x ) dx − xS(x)f (t, x) x

= xQbreak (t, x), where we have used Eq. (4) and Eq. (5).

(A.1) x 0

 y b(t, y , x) dy

= x according to

References Abramowitz, M., Stegun, I.A., 1972. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover, New York. Daniele, L.M., Dennis, R.V., Rodney, O.F., 2003. Quadrature method of moments for aggregation–breakage processes. J. Colloid Interface Sci. 258, 322–334. Everson, R.C., Eyre, D., Campbell, Q.P., 1997. Spline method for solving continuous batch grinding and similarity equations. Comput. Chem. Eng. 21, 1433–1440.

2095

Filbet, F., Laurençot, P., 2004. Numerical simulation of the Smoluchowski coagulation equation. SIAM J. Sci. Comput. 25 (6), 2004–2048. Gunawan, R., Fusman, I., Braatz, R.D., 2004. High resolution algorithms for multidimensional population balance equations. A.I.Ch.E. J. 50, 2738–2749. Gunawan, R., Fusman, I., Braatz, R.D., 2008. Parallel high-resolution finite volume simulation of particulate processes. A.I.Ch.E. J. 54, 1449–1458. Hill, P.J., Ng, K.M., 1995. New discretization procedure for the breakage equation. A.I.Ch.E. J. 41, 11204–11216. Hounslow, M.J., Pearson, J.M.K., Instone, T., 2001. Tracer studies of high shear granulation: population balance modeling. A.I.Ch.E. J. 47, 1984–1999. Hulbert, H.M., Katz, S., 1964. Some problems in particle technology. Chem. Eng. Sci. 19, 555–574. Koren, B., 1993. A robust upwind discretization method for advection, diffusion and source terms. In: Vreugdenhill, C.B., Koren, B. (Eds.), Numerical Methods for Advection–Diffusion Problems, vol. 45. Vieweg Verlag, Braunschweig, pp. 117–138 (Notes on Numerical Fluid Mechanics, chapter 5). Kumar, J., Numerical approximations of population balance equations in particulate systems. Ph.D. Thesis, Faculty of Mathematics, Otto-von-Guericke University, 2006. ¨ L., 2006. Improved accuracy Kumar, J., Peglow, M., Warnecke, G., Heinrich, S., Morl, and convergence of discretized population balances: the cell average technique. Chem. Eng. Sci. 61, 3327–3342. Kostoglou, M., Karabelas, A.J., 2002. An assessment of low-order methods for solving the breakage equation. Powder Technol. 127, 116–127. Kostoglou, M., Karabelas, A.J., 2004. Optimal low order methods of moments for solving the fragmentation equation. Powder Technol. 143–144, 116–127. Kumar, S., Ramkrishna, D., 1996. On the solution of population balance equations by discretization—I. A fixed pivot technique. Chem. Eng. Sci. 51, 1311–1332. Kumar, S., Ramkrishna, D., 1997. On the solution of population balance equations by discretization—III. Nucleation, growth and aggregation of particles. Chem. Eng. Sci. 52, 4659–4679. Lage, P.L.C., 2002. Comments on the “An analytical solution to the population balance equation with coalescence and breakage—the special case with constant number of particles”, by D.P. Patil and J.R.G Andrews. Chem. Eng. Sci. 57, 4253–4254. Lee, K., Matsoukas, T., 2000. Simultaneous coagulation and breakage using constant—N Monte Carlo. Powder Technol. 110, 82–89. LeVeque, R.J., 2002. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, New York, NY. 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. Chem. Eng. Sci. 57, 3715–3732. Ma, A., Tafti, D.K., Braatz, R.D., 2002. High resolution simulation of multidimensional crystal growth. Ind. Eng. Chem. Res. 41, 6217–6223. Mishra, B.K., 2000. Monte Carlo simulation of particle break process during grinding. Powder Technol. 110, 246–252. Nicmanis, M., Hounslow, M.J., 1998. Finite element methods for steady-state population balance equations. A.I.Ch.E. J. 44, 2258–2272. Perry, R.-H., Green, D.W., 1997. Perry's Chemical Engineers' Handbook. 7th ed. McGraw-Hill, New York. 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. Comput. Chem. Eng. 30, 1119–1131. Qamar, S., Ashfaq, A., Elsner, M.P., Angelov, I., Warnecke, G., Seidel-Morgenstern, A., 2007. Adaptive high resolution schemes for multidimensional population balances in crystallization processes. Comput. Chem. Eng. 31, 1296–1311. Qamar, S., Warnecke, G., 2006. Solving population balance equations for twocomponent aggregation by a finite volume scheme. Chem. Eng. Sci. 62, 679–693. Qamar, S., Warnecke, G., 2007. Numerical solution of population balance equations for nucleation, growth and aggregation processes. Comput. Chem. Eng. 31, 1576–1589. Ramkrishna, D., 1985. The status of population balances. Rev. Chem. Eng. 3 (1), 49–95. Ramkrishna, D., 2000. Population Balances: Theory and Applications to Particulate Systems in Engineering. Academic Press, San Diego, CA. Randolph, A., Larson, M.A., 1988. Theory of Particulate Processes. 2nd ed. Academic Press, San Diego, CA. Vanni, M., 1999. Discretization procedure for the breakage equation. A.I.Ch.E. J. 45, 916–919. Woo X.Y., Tan, R.B.H., Chow, P.S., Braatz, R.D., 2006. Simulation of mixing effects in antisolvent crystallization using a coupled CFD–PDF–PBE approach. Crystal Growth Des. 1291–1303. Ziff, R.M., McGray, E.D., 1985. The kinetics of cluster fragmentation and depolymerization. J. Phys. A Math. Gen. 18, 3027–3037.