Lattice Boltzmann method for population balance equations with simultaneous growth, nucleation, aggregation and breakage

Chemical Engineering Science 69 (2012) 316–328

Contents lists available at SciVerse ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

Lattice Boltzmann method for population balance equations with simultaneous growth, nucleation, aggregation and breakage Aniruddha Majumder a, Vinay Kariwala a,n, Santosh Ansumali b, Arvind Rajendran a a b

School of Chemical & Biomedical Engineering, Nanyang Technological University, Singapore 637459, Singapore Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur, Bangalore 560 064, India

a r t i c l e i n f o

a b s t r a c t

Article history: Received 4 June 2011 Received in revised form 4 October 2011 Accepted 21 October 2011 Available online 4 November 2011

Lattice Boltzmann method (LBM) is developed for solution of one-dimensional population balance equations (PBEs) with simultaneous growth, nucleation, aggregation and breakage. Aggregation and breakage, which act as source terms in PBEs, are included as force terms in LBM formulation. The force terms representing aggregation and breakage are evaluated by fixed pivot (FP) method. Multiscale analysis is used to derive the kinetic equations associated with LBM, whose long-time large-scale solution provides the solution of the PBE. A coordinate transformation is proposed, which allows the use of non-uniform grid for LBM to obtain accurate solution of PBE with moderate number of grid points. The performance of the proposed LBM-FP method is compared with finite volume (FV) and method of characteristics (MOC) combined with FP (MOC-FP) methods. Using benchmark examples, the proposed LBM-FP method is shown to be useful for solving PBEs due to its computational efficiency and ability to handle a wide range of problems. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Aggregation Breakage Dynamic simulation Lattice Boltzmann method Particulate process Population balance

1. Introduction Population balance modeling has become an important tool to model particulate processes in various disciplines of science and engineering, e.g., crystallization (Randolph and Larson, 1962; Ramkrishna, 2000), granulation (Poon et al., 2009), milling (Bilgili and Scarlett, 2005), emulsion polymerization (Immanuel et al., 2002), aerosols (Ramabhadran et al., 1976) and biological systems (Mantzaris et al., 2001). As the name signifies, population balance seeks to describe the behavior of the population of the particles by tracking the evolution of the number density function defined in the particle state space. Population balance equation (PBE) takes into account various mechanisms, such as growth, nucleation, aggregation and breakage, by which particles of a particular state can either form or disappear from the system. In this work, we assume that the particles are uniformly distributed in the well-mixed control volume, which implies that the particle size distribution (PSD) is independent of the spatial coordinate. Moreover, only one characteristic property of the particles is considered. The resulting one-dimensional (1D) PBE takes the following form: @ @ nðx,tÞ þ ðGðx,tÞnðx,tÞÞ ¼ Qnuc þ Qagg þ Qbreak , @t @x n

Corresponding author. Tel.: þ65 6316 8746; fax: þ65 6794 7553. E-mail address: [email protected] (V. Kariwala).

0009-2509/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2011.10.051

ð1Þ

where nðx,tÞ is the PSD, Gðx,tÞ represents the growth rate of the particle at state x and the terms on the right hand side of Eq. (1) denote the rates of nucleation (nuc), aggregation (agg) and breakage (break) given as Qnuc ¼ B0 ðtÞdðx0 Þ, Qagg ¼

ð2Þ

Z 1 Z 1 x 0 0 nðxx0 ,tÞnðx0 ,tÞaðxx0 ,x0 Þ dx  nðx,tÞnðx0 ,tÞaðx,x0 Þ dx , 2 0 0 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Birth term due to aggregation, RBa

Qbreak ¼

Z

1

0

Death term due to aggregation, RDa

bðx,x0 ÞGðx0 Þnðx0 ,tÞ dx  x |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Death Birth term due to breakage RBb

ð3Þ

GðxÞnðx,tÞ |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl}

,

ð4Þ

term due to breakage RDb

where B0 is the nucleation rate, d is the Dirac delta function, x0 is the size of the nuclei, aðx,x0 Þ is the aggregation frequency denoting the rate at which two particles of size x and x0 are combined together to produce a particle of size x þ x0 , GðxÞ is the breakage rate and bðx,x0 Þ denotes the probability of formation of particle of size x due to breakage of particle of size x0 . In general, the PBE in Eq. (1) has to be supplemented with appropriate equations for the kinetics (e.g., growth and nucleation) as well as mass and energy balances (Lang et al., 1999; Myerson, 2002). PBEs are hyperbolic partial differential equations. In most practical cases, they do not have analytical solution and hence need to be solved numerically. During the past few decades, a number of methods have been developed for numerical solution

A. Majumder et al. / Chemical Engineering Science 69 (2012) 316–328

of PBEs. Several reviews of these techniques are available in the literature (Ramkrishna, 2000; Motz et al., 2002; Nopens et al., 2005; Costa et al., 2007; Kumar et al., 2009; Silva et al., 2010). In the subsequent discussion, we provide an overview of the relevant methods for solving PBEs involving simultaneous growth, nucleation, aggregation and breakage. Some of the available methods are used to track the evolution of the moments of the PSD. Various quadrature method of moments (QMOM) are used for this purpose (Marchisio et al., 2003; Marchisio and Fox, 2005; Gimbun et al., 2009; Qamar et al., 2011; Kariwala et al., in press). However, full PSD is often required for design and control purposes (Nagy, 2009) and reconstruction of the PSD from the moments can be difficult (Aamir et al., 2009). Finite element method (FEM) has been used for solving PBEs with simultaneous growth and aggregation to obtain information about the full PSD (Nicmanis and Hounslow, 1998; Mahoney and Ramkrishna, 2002; Alexopoulos et al., 2009). Although accurate, FEM is usually computationally more expensive than other methods, e.g., discretized methods. Monte Carlo (MC) techniques are also used for solving PBEs with simultaneous processes (van Peborgh Gooch and Hounslow, 1996; Smith and Matsoukas, 1998; Zhao et al., 2007). The major drawbacks of this method are that the accuracy largely depends on the number of sample particles considered (Meimaroglou and Kiparissides, 2007) and the involved computational expense is relatively high for low dimensional PBEs. In recent years, another class of methods, known as discretized methods, have become popular for solving PBEs due to their efficiency and accuracy. These methods include fixed pivot (FP) method (Kumar and Ramkrishna, 1996a, 1997; Lim et al., 2002), moving pivot (MP) method (Kumar and Ramkrishna, 1996b), cell average technique (Kumar et al., 2006, 2008b), hierarchical two-tier method (Immanuel and Doyle, 2003, 2005), two-level discretization algorithm (Pinto et al., 2007, 2008) and finite volume (FV) method (Qamar and Warnecke, 2007; Qamar et al., 2009). The FP method, proposed by Kumar and Ramkrishna (1996a), has been widely used for solving PBEs with aggregation and breakage due to its easy implementation and ability to handle a variety of grids. In this method, the computational domain is discretized into small cells. The size of all the particles in each cell is represented by a size that falls within that cell and is called the pivotal point. The key idea in this method is to distribute the newly born particles, which do not fall on the pivotal points, among the neighboring pivotal points such that some properties of interest are preserved exactly. When a coarse grid is used, the FP method overpredicts the PSD for aggregation process in the region where the PSD changes rapidly. In order to overcome this problem, MP method was proposed (Kumar and Ramkrishna, 1996b). In this method, the pivotal size in each cell is adjusted according to the non-uniformity of the PSD in that cell (Kumar and Ramkrishna, 1996b). However, the MP method is more complex than the FP method and the resulting ODEs are difficult to solve (Kumar et al., 2006). Since these methods are developed for pure aggregation and breakage processes, they need to be combined with other methods for handling simultaneous processes, such as growth and aggregation. Such methods include combination of method of characteristics (MOC) and FP method (Kumar and Ramkrishna, 1997; Lim et al., 2002) and combination of weighted essentially non-oscillatory (WENO) and FP method (Lim et al., 2002). The cell average technique, developed initially for pure aggregation and breakage processes, is similar to the FP method but with the improved assignment rule for the newly born particles to the neighboring cells (Kumar et al., 2006, 2008a). This method has been extended to solve simultaneous processes including growth by reformulating the growth term (Kumar et al., 2008b). The two-tier technique comprises of two steps. The first step involves the calculation of the rates of nucleation, growth and aggregation. This information is then used in the second step to

317

update the PSD (Immanuel and Doyle, 2003, 2005). A variant of the two-tier technique is the two-level discretization technique (Pinto et al., 2007, 2008). In this technique, a coarse mesh is used for aggregation, while a fine mesh is used for growth and nucleation based on the observation that the effect of aggregation on PSD is less sensitive to mesh size than the effects of nucleation and growth. In FV method, the PBE is reformulated in a mass conservation form so that the aggregation and breakage terms can be written as flux terms (Filbet and Laurencot, 2004; Qamar and Warnecke, 2007; Qamar et al., 2009). The reformulated PBE can be easily solved using FV method available for conservation laws (LeVeque, 2002). Many methods mentioned earlier are not very efficient for model based online optimization and control, where repeated solution of the model equations is required. Efficient solution of the PBE is also crucial when computational fluid dynamics (CFD) is coupled with PBE to take into account the effect of non-ideal mixing (Marchisio and Fox, 2005; Woo et al., 2009). Recently, we have proposed a mesoscopic particle based method, namely lattice Boltzmann method (LBM), for efficiently solving PBEs with growth and nucleation (Majumder et al., 2010a, in press). In this paper, LBM is extended for efficient solution of 1D PBE with simultaneous growth, nucleation, aggregation and breakage phenomena. This extension requires incorporation of aggregation and breakage processes as force terms in the LB formulation. This force term is evaluated by applying the FP method and thus the resulting method is referred to as LBM-FP method in the subsequent discussion. Detailed multiscale analysis is carried out to derive the kinetic equations, whose long-time large-scale solution provides the solution of the desired PBE. The traditional LBM is only applicable with uniform grid. The use of such grid requires a large number of grid points to capture the PSD with reasonable resolution for aggregation and breakage processes and hence make the numerical scheme computationally expensive. In order to overcome this problem in the case of FP method, the use of non-uniform grid (e.g., geometric grid) has been proposed by Kumar and Ramkrishna (1996a). However, the implementation of such non-uniform grid in LBM requires interpolation of the discrete Boltzmann distribution function (He et al., 1996). In this paper, a coordinate transformation is proposed, which allows the use of non-uniform grid (e.g., logarithmic grid) without performing any interpolation and thus avoids the inaccuracy and computational cost resulting due to the interpolation. The performance of the proposed LBM-FP technique is assessed for pure aggregation and breakage processes, as well as simultaneous processes, such as processes with simultaneous growth, nucleation, aggregation and breakage phenomena. The accuracy and efficiency of the proposed scheme is compared with the FP (Kumar and Ramkrishna, 1996a), MOC-FP (Kumar and Ramkrishna, 1997) and FV (Qamar et al., 2009) methods. It is found that the LBM-FP technique has similar accuracy as FP and FV methods for pure aggregation and pure breakage. However, the proposed technique is found to be more accurate and more efficient than the FV method for the growth dominant simultaneous processes, e.g., simultaneous growth and aggregation. In comparison to MOC-FP method, the proposed method is found to be more accurate when nucleation is present and when growth rate is constant. The rest of the article is organized as follows: the principle of LBM is discussed and detailed derivation of LBM for solving advection equation with source term is presented in Section 2. In Section 3, adaptation of the LBM for solving PBEs with simultaneous growth, nucleation, aggregation and breakage is discussed by identifying the analogy between the advection equation and the PBE. Subsequently, some case studies are presented in Section 4, where the performance of the proposed LBM-FP method is compared with FP, MOC-FP and FV methods. Finally, some concluding remarks are presented in Section 5.

318

A. Majumder et al. / Chemical Engineering Science 69 (2012) 316–328

2. Lattice Boltzmann method

Right moving

Lattice based methods (Frisch et al., 1986; McNamara and Zanetti, 1988; Higuera et al., 1989) have received significant attention in the past two decades for simulation of hydrodynamics. The major difference between lattice based and other traditional methods, e.g., finite difference method, is that most of the traditional methods follow a top-down approach, where the governing partial differential equations (PDEs) at the macroscopic scale are discretized to obtain ordinary differential equations (ODEs) or algebraic equations. On the other hand, lattice based methods adopt a bottom-up approach, where instead of using discretization, pseudo particle kinetics is used in such a way that the governing equations are recovered at the appropriate scale. LBM was initially developed as a spin-off of one such method known as Lattice Gas Cellular Automata (LGCA) (Frisch et al., 1986; Succi, 2001). The key distinguishing feature of LBM is the introduction of Boltzmann equation to overcome the problem of statistical noise encountered in LGCA. With its roots in kinetic theory, LBM uses a simplified representation of the microscopic state of the process so that the model is computationally viable and can capture the essential description of the system at the macroscopic level. LBM has now established itself as a popular scheme due to its ability to provide fast, easily implementable code and versatility (Aidun and Clausen, 2009). It has been applied to various problems including hydrodynamics (Chen and Doolen, 1998), turbulence (Chen et al., 2003), multi-phase flow (Mantle et al., 2001), microflow (Ansumali et al., 2007; Kim et al., 2008; Yudistiawan et al., 2010), non-Newtonian fluids (Singh et al., 2011) and crystallization (Majumder et al., 2010a, in press). Several good reviews of LBM are available in the literature (Benzi et al., 1992; Chen and Doolen, 1998; Succi, 2001; Aidun and Clausen, 2009). PBE with aggregation and breakage is analogous to the advection equation with source term. Thus, in the subsequent discussion, we first derive LBM for the advection equation with source term and then apply the developed scheme for solving PBEs with aggregation and breakage by identifying their similarities later in the paper. 2.1. Principle The advection equation with a source term can be written as @r @ðvrÞ ¼ F, þ @r @t

ð5Þ

where r is the concentration of the species being transported (passive scalar), v is the advection velocity and F is the source term. In order to solve this advection equation, some fictitious particles, which resemble groups of molecules, are considered at the mesoscopic scale. While molecules move randomly in space, these fictitious particles are restricted to move with certain velocities, which are carefully chosen such that they are consistent with the symmetry and isotropy requirement of macroscopic dynamics (Succi, 2001). Finding such a set of velocities involves a trial and error procedure. One typically starts with a set of velocities on a lattice based on computational considerations, and checks whether it is possible to construct a model, whose long-term large-scale solution matches with the macroscopic equation (Succi, 2001). In the discrete form, the simplest possible model has three types of fictitious particles with velocities ci ¼ f0,c,cg, i.e., stationary (0), right moving (þ), left moving () (Karlin et al., 2006), as illustrated in Fig. 1. In order to solve Eq. (5) with LBM, the following kinetic equation with force (source) term is used (Dawson et al., 1993): @f i @f 1 eq þ ci i ¼  ðf i f i Þ þ F i , @t @r t

i ¼ 0; 1,2,

ð6Þ

Stationary Left moving Fig. 1. 1D lattice model showing fictitious particles with three velocities.

where t 4 0 is the relaxation time, fi is the discrete Boltzmann eq distribution function, f i is the equilibrium discrete Boltzmann distribution function, ci is the discrete velocity and Fi is the discrete force term. In the simplest form, Fi can be taken as Fi ¼wiF (Dawson et al., 1993), where wi is the weight associated with velocity ci. In Eq. (6), the subscripts 0, 1 and 2 refer to stationary, right moving and left moving particles, respectively. The left hand side of Eq. (6) denotes the free flight, while the right hand side represents the relaxation of the particles to equilibrium (Bhatnagar–Gross–Krook (BGK) approximation (Bhatnagar et al., 1954) for collision). The equilibrium distributions can be found by minimizing appropriate entropy function subject to the constraints of mass and momentum conservation. The discrete form of the relevant entropy function, also known as H function, is given as (Karlin et al., 1999; Ansumali and Karlin, 2002a)   2 X f H¼ f i ln i 1 : ð7Þ wi i¼0 For the 1D advection equation, the weights can be selected as w0 ¼ 4=6, w1 ¼ 1=6, w2 ¼ 1=6 (Karlin et al., 1999; Ansumali and Karlin, 2002b). The constraints for local conservation of mass and momentum are given as 2 X

f i ¼ f 0 þ f 1 þ f 2 ¼ r,

ð8Þ

ci f i ¼ cðf 1 f 2 Þ ¼ ru,

ð9Þ

i¼0 2 X i¼0

where u is the local average velocity. By solving the minimization problem, we get the equilibrium values of the discrete Boltzmann eq distribution function f i as (Ansumali and Karlin, 2000)  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2r eq 2 1þ u2 =c2s , ð10Þ f 0 ðr,uÞ ¼ 3 eq

f 1 ðr,uÞ ¼ eq

f 2 ðr,uÞ ¼

r 3

r 3

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðucc2s Þ=2c2s þ 1þ u2 =c2s ,

ð11Þ

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðuc þc2s Þ=2c2s þ 1 þu2 =c2s ,

ð12Þ

pffiffiffi where cs ¼ c= 3 is analogous to the speed of sound in the system. Conserving both mass and momentum during collision will recover the Navier–Stokes equation in the long-time long-scale limit, whereas conserving only mass will recover the advection equation (Karlin et al., 2006; Majumder et al., 2010a). Since we are interested in the advection equation, equilibrium distribution eq f i ðrðf Þ,uÞ is calculated at the given velocity u ¼v (Karlin et al., 2006). It is shown next that breaking this momentum conservation leads to advection–diffusion equation for r in the long-time and large-scale limit. 2.2. Macroscopic equation from kinetic equation with source term In this section, multiscale analysis, known as Chapman–Enskog analysis, is used to verify whether the kinetic equation in Eq. (6) is consistent with the advection equation at the macroscopic level.

A. Majumder et al. / Chemical Engineering Science 69 (2012) 316–328

Using Eqs. (8) and (9) the kinetic equation given by Eq. (6) can be equivalently written in terms of moment chain as 0th moment :

@r @ðruÞ ¼ F, þ @r @t

ð13Þ

1st moment :

@ðruÞ @P 1 þ ¼  ðrurvÞ, @t @r t

ð14Þ

2nd moment :

@P @ðruÞ 1 1 þ c2 ¼  ðPPeq Þ þ c2 F, @t @r 3 t

ð15Þ

where P¼

2 X

c2i f i ,

ð16Þ

i¼0

P eq ¼

2 X

eq

c2i f i :

ð17Þ

i¼0

In order to obtain the macroscopic PDE from the LBM formulation, it is necessary to formally separate different time scales. For this purpose, the non-conserved variables, i.e., u and P, are expanded around their equilibrium values in terms of t. The time derivative operator is also expanded in terms of t as follows: ð1Þ

u ¼ v þ tu

2 ð2Þ

þt u

þ ,

ð18Þ

P ¼ P eq þ tPð1Þ þ t2 Pð2Þ þ    ,

ð19Þ

@ @ð0Þ @ð1Þ ¼ þt þ , @t @t @t

ð20Þ

where advection velocity v is the equilibrium value of the local average velocity u and the superscripts within the parenthesis denote the order of the expansion. Upon substituting these expansions in Eq. (13), we have  ð0Þ  @ @ð1Þ @ þt þ    r þ ½rv þ truð1Þ þ    ¼ F: ð21Þ @r @t @t Taking Oð1Þ terms @ð0Þ @ r ¼  ðrvÞ þ F: @r @t

ð22Þ

t

ð23Þ

The Oð1Þ terms of Eq. (23) are ruð1Þ ¼

@ð0Þ @ ðrvÞ þ Peq : @r @t

ð24Þ

By substituting the expression for Peq and using Eq. (22), Eq. (24) can be written as    @ð0Þ @ v2 ruð1Þ ¼ v r þ c2s r 1þ 2 , ð25Þ @r @t cs ruð1Þ ¼ v ruð1Þ ¼ rv

@ @ @v ðrvÞ þ vF þ ðc2s þ v2 Þ r þ2rv , @r @r @r

ð26Þ

@v @r þ c2s þvF: @x @r

ð27Þ

ð28Þ

ð29Þ

In Eq. (30), D ¼ tc2s is the diffusion coefficient, which can be made small by choosing sufficiently small relaxation time t. The last two terms on the right hand side of Eq. (30) are the anomaly terms. The last term arises due to the space dependence of the advection velocity. When v is a smooth function of space coordinate r, this anomaly term is OðtMa2 Þ, where Ma ¼ v=cs is the Mach number. The contribution from this term can be neglected for small Ma, e.g., Ma  0:01, and vanishes when v is independent of r. Note that Ma cannot be chosen to be very small, as smaller Ma (i.e., larger cs) increases the diffusion coefficient D. The anomaly term involving the force ðtð@vF=@rÞÞ can be nullified by considering the following correction term in the kinetic equation: @f i @f 1 eq þ ci i ¼  ðf i f i Þ þwi F þ wi ci C, @t @r t

ð31Þ

where C is the unknown parameter included in the correction term. Now, following the same steps for the Chapman–Enskog analysis as used earlier, the following equation is recovered at the macroscopic level:   @ @ @2 r @vF @C @ @v tc2s : ð32Þ þt r þ ðrvÞ ¼ F þD 2 þ t rv @t @r @r @r @r @r @r In Eq. (32), the third term on the right hand side is the anomaly term involving F and the fourth term is the correction term. Thus, the condition to nullify this anomaly term is @ ½tvFtc2s C ¼ 0, @r

ð33Þ

) tvF ¼ tc2s C þ C,

ð34Þ

v F, c2s

ð35Þ

where the integration constant C has been taken as zero. With this choice, the kinetic equation with the correction term becomes   @f i @f 1 cv eq þ ci i ¼  ðf i f i Þ þwi F 1 þ i2 : ð36Þ @t @r t cs Effectively, to nullify the anomaly term due to F in Eq. (30), instead of choosing the force term in the kinetic equation as wiF as suggested by Dawson et al. (1993), one needs to use   cv ð37Þ F i ¼ wi F 1 þ i2 : cs 2.3. Discretization of the kinetic equation In order to solve the kinetic equation numerically, we integrate it over Dt using the trapezoidal rule (which is OðDt 2 Þ accurate) to obtain f i ðr þ ci Dt,t þ DtÞ  f i ðr,tÞ

Then, the expression for momentum ru in linear order of t becomes

ru ¼ rv þ truð1Þ þ    ,

@v @r tc2s tvF: @r @r

By substituting the expression for ru in Eq. (13), we obtain   @ @ @2 r @vF @ @v þt : ð30Þ r þ ðrvÞ ¼ F þD 2 þ t rv @t @r @r @r @r @r

) C¼

Similarly, by substituting the expansion values in Eq. (14), the following expression is obtained:  ð0Þ  @ @ð1Þ @ þt þ    ½rv þ truð1Þ þ    þ ½P eq þ tPð1Þ þ    @r @t @t 1 ¼  ½rv þ truð1Þ þ    rv:

ru ¼ rvtrv

319

Dt

eq

½ðf ðr,tÞf i ðr,tÞÞ 2t i eq þ ðf i ðr þ ci Dt,t þ DtÞf i ðr þ ci Dt,t þ DtÞÞ

þ

Dt 2

½F i ðr þ ci Dt,t þ DtÞ þ F i ðr,tÞ:

ð38Þ

It can be noted that the scheme given by Eq. (38) is not explicit as f i ðr þci Dt,t þ DtÞ appears on both sides. In order to obtain an

320

A. Majumder et al. / Chemical Engineering Science 69 (2012) 316–328

explicit numerical scheme, we define a transformation as f^ i ¼ f i þ

Dt 2t

eq

ðf i f i Þ

Dt 2

Fi:

ð39Þ

With this transformation, Eq. (38) can be rewritten in terms of f^ i as

Dt eq f^ i ðr þ ci Dt,t þ DtÞ  f^ i ðr,tÞ ðf i ðr,tÞf i ðr,tÞÞ þ DtF i ðr,tÞ:

t

ð40Þ

Based on Eqs. (8) and (39) the density r in terms of f^ i can be found as 2 X

Dt f^ i ¼ r 2 i¼0 ) r¼

X Fi,

ð41Þ

i

2 X

Dt F: f^ i þ 2 i¼0

ð42Þ eq

Moreover, from Eqs. (10)–(12), f i is a functional of fi and depends on the particle population only via its lower order moments r and ru. Substituting for r in Eqs. (10)–(12) with the expression in Eq. (42), it is found that at the same local average velocity u, the eq eq equilibrium distributions of f i and f^ i are related as eq Dt eq F: f^ i ¼ f i  2 i

ð43Þ

Now, from Eqs. (39) and (43) eq

f i f i ¼

eq 2t ðf^ f^ Þ: 2t þ Dt i i

ð44Þ

eq

Substituting for ðf i f i Þ from Eq. (44) in Eq. (40), we obtain   2Dt 2Dt ^ eq f^ ðr,tÞ þ f ðr,tÞ þ DtF i , f^ i ðr þ ci Dt,t þ DtÞ  1 2t þ Dt i 2t þ Dt i ð45Þ eq ) f^ i ðr þci Dt,t þ DtÞ  ð1abÞf^ i ðr,tÞ þ abf^ i ðr,tÞ þ DtF i ,

ð46Þ

where a ¼ 2 and b ¼ Dt=ð2t þ DtÞ. Thus, we arrive at an explicit expression in terms of f^ i . Although f^ i are not exactly particle populations, they serve the purpose of solving the advection equation since the quantity r can be recovered from f^ i using Eq. (42). Numerical implementation of the LBM scheme developed for advection equation requires solving Eq. (46). In practice, this equation is split into two consecutive steps: (1) collision; and (2) streaming. The collision step involves evaluation of the right hand side of Eq. (46), while streaming step involves shifting the values of the distribution functions to the next lattice node based on their discrete velocities.

and Ramkrishna, 1996a). On the other hand, a non-uniform grid such as geometric or logarithmic grid can provide better resolution of the size distribution by covering the entire size range of the particles with moderate number of grid points (Batterham et al., 1981; Kumar and Ramkrishna, 1996a). However, in such a grid, the newly born particles due to aggregation and breakage may have sizes which do not fall on the grid points. This problem can be overcome by redistributing the newly born particles among the neighboring grid points so that some properties of interest are preserved. We use the fixed pivot (FP) method proposed by Kumar and Ramkrishna (1996a) for this purpose, where the zeroth and the first order moments of the PSD, i.e., total number and total size of the particles are considered to be conserved when particle size is the independent variable. The LBM discussed earlier for solving the advection equation requires a uniform grid for implementation, as the grid spacing should be such that during the streaming step all the particles from one node reach the next node. Non-uniform grid can be implemented with an interpolation based LBM, where the populations that do not reach the next node in the streaming step are interpolated (He et al., 1996; Chikatamarla and Babu, 2005). In order to avoid the inaccuracy and computational cost due to interpolation, we propose a coordinate transformation, which usually leads to a logarithmic grid in the original coordinate system (Majumder et al., 2010b). Similar to geometric grid, the logarithmic grid allows capturing the PSD with good resolution with moderate number of grid points. The proposed coordinate transformation is defined as z ¼ zðxÞ. As the number of particles in a specified size interval in both the coordinate systems are the same, the PSDs in terms of z and x are related as (Roussas, 2003)   dx hðz,tÞ ¼ nðx,tÞ , ð47Þ dz where hðz,tÞ is the number density in the transformed coordinate system. We consider the case of size-dependent growth rate of the particles, which is usually modeled as (Abegg et al., 1968; Lim et al., 2002; Gunawan et al., 2004) Gðx,tÞ ¼ FðtÞOðxÞ,

where FðtÞ and OðxÞ characterize the size and time dependence of the growth rate. We assume that OðxÞ is continuous and positive over the domain ½xmin ,xmax  to make sure that the mapping from one coordinate system to another is unique and to avoid singularity (Majumder et al., 2010b). Substituting for nðx,tÞ and Gðx,tÞ on the left hand side of Eq. (1), we obtain  1 !  1 ! @ dx @ dx h OðxÞh þ FðtÞ ¼ Qnuc þ Qagg þ Qbreak : @t dz @x dz

3. Adaptation of LBM for solving PBEs The PBE given by Eq. (1) is analogous to the advection equation given by Eq. (5), where the PSD n is analogous to the fluid density r, the growth rate G is analogous to the local fluid velocity v, the space coordinate r is analogous to the characteristic size x and the terms on the right hand side of the PBE (Qnuc , Qagg and Qbreak ) are analogous to the source term F. Thus, the LBM scheme discussed in the previous section for advection equation with source term F can be applied to solve PBEs with simultaneous growth, nucleation, aggregation and breakage. While solving such a PBE, the nucleation term is implemented as a boundary condition as the size of the nuclei is assumed to be represented by the smallest grid considered. The integrals related to aggregation and breakage can be evaluated using numerical integration techniques easily, if a uniform grid is used. However, a uniform grid requires very large number of grid points in order to cover the entire size range of the particles with reasonable accuracy and is thus computationally expensive (Kumar

ð48Þ

ð49Þ Without loss of generality, we may define dx ¼ O ðxÞ, dz

ð50Þ

) z ¼ cðxÞ,

ð51Þ

where O ðxÞ is a function of x that defines the coordinate transformation and Z dx þ k, ð52Þ c¼ O ðxÞ where k is the constant of integration. Thus, in the transformed coordinate z, Eq. (1) becomes !!  1 1 dx @h @ Oðc ðzÞÞ þ FðtÞ h ¼ Qnuc þQagg þ Qbreak : ð53Þ dz @t @z O ðc1 ðzÞÞ

A. Majumder et al. / Chemical Engineering Science 69 (2012) 316–328

Now, substituting for nðx,tÞ and x on the right hand side of Eq. (53), we get 1

Qnuc ¼ B0 dðc

ðz0 ÞÞ:

ð54Þ

Similarly, other terms on the right hand side of Eq. (53) become  1 Z z 1 dx 1 1 1 0 aðc ðzÞc ðz0 Þ, c ðz0 ÞÞhðzz0 ,tÞhðz0 ,tÞ dz , RBa ¼ 2 dz 0 ð55Þ RDa ¼



1 Z 1 dx 1 1 0 aðc ðzÞ, c ðz0 ÞÞhðz,tÞhðz0 ,tÞ dz , dz 0

1 Z 1 dx 1 1 1 0 bðc ðzÞ, c ðz0 ÞÞGðc ðz0 ÞÞhðz0 ,tÞ dz , RBb ¼ dz z 

1 dx Gðc1 ðzÞÞhðz,tÞ: dz

ð57Þ

ð58Þ

Q nuc

Z 1 z 1 1 1 0 þ aðc ðzÞc ðz0 Þ, c ðz0 ÞÞhðzz0 ,tÞhðz0 ,tÞ dz 2 0 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

9nei ni 9Dxi ,

ð65Þ

ð66Þ

i¼1

1

1

0

1

bðc ðzÞ, c ðz0 ÞÞGðc ðz0 ÞÞhðz0 ,tÞ dz Gðc ðzÞÞhðz,tÞ : |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} z |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} R Db

R Bb

ð59Þ For processes with size-dependent growth, we can choose

O ðxÞ ¼ OðxÞ. In this case, we usually have a logarithmic grid in the original coordinate system, since OðxÞ is often modeled as OðxÞ ¼ 1 þ gx, where g 4 0 is a parameter that determines the extent of size dependence (Gunawan et al., 2004). With this transformation, Eq. (59) becomes a PBE with size-independent growth rate, i.e., @h @h þ FðtÞ ¼ Q nuc þ R Ba þ R Da þ R Bb þ R Db : @t @z

ð60Þ

The developed LBM scheme now can be used efficiently for solving Eq. (60). On the other hand, if the growth rate is size independent, e.g., OðxÞ ¼ 1, we can define O ðxÞ ¼ 1 þ gx. In such case, Eq. (59) simplifies as @h @ þ ðGhÞ ¼ Q nuc þ R Ba þ R Da þ R Bb þR Db , @t @z

ð61Þ

where Gðz,tÞ ¼

N X

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N uX ðnei ni Þ2 Dxi , L2 -norm ¼ t

R Da

1

In this section, numerical examples are considered for pure breakage, pure aggregation, simultaneous growth–aggregation and simultaneous growth–nucleation–aggregation–breakage processes. The problems are solved using the proposed LBM-FP method and the results are compared with FP (Kumar and Ramkrishna, 1996a), MOC-FP (Kumar and Ramkrishna, 1997) and FV methods (Qamar et al., 2009) in terms of error norms and computation time. The L1 - and L2 -norms of the error are defined as follows:

i¼1

R Ba

1

ð64Þ

where f^ i,buff is the population density in the buffer point considered.

L1 -norm ¼

Z 1 1 1 0  aðc ðzÞ, c ðz0 ÞÞhðz,tÞhðz0 ,tÞ dz 0 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} þ

eq f^ i,buff ¼ f^ i ðhðz0 ,tÞÞ,

4. Numerical examples

Finally, by combining Eqs. (53)–(58), we obtain ! 1 @h @ Oðc ðzÞÞ 1 1 þ FðtÞ h ¼ B0 dðc ðz0 ÞÞO ðc ðzÞÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} @t @z O ðc1 ðzÞÞ

Z

condition can be implemented by considering a buffer point before the first grid point in the computational domain and continually refilling this point with the equilibrium populations corresponding to the number density at size z0 as (Succi, 2001; Majumder et al., 2010a)

ð56Þ



RDb ¼

321

FðtÞ : expðgzÞ

where xi are the discrete points of the chosen grid, Dxi is width of the ith cell and nei is the exact solution. For simulation of pure aggregation and breakage processes, the LBM parameter c is chosen such that the time step Dt ¼ Dx=c  0:01 s. Selecting the relaxation parameter b closer to 1 provides smaller diffusion coefficient but it may introduce oscillations, especially for nonsmooth distributions. Thus b is chosen such that the diffusion coefficient is sufficiently small and oscillations are minimal. A typical value of b used here is 0.99. For FV, MOC-FP and pure FP methods, semidiscrete formulation is used and the Matlab ODE solver ode45 is applied to solve the resulting ODEs. The default values of absolute tolerance of 106 and relative tolerance of 103 are used for the ODE solver unless otherwise stated. All computations are performed using a Windows Vista PC with an Intel CoreTM 2 Duo s Processor 6600 (2.40 GHz, 4 GB RAM) using Matlab 2007a.

4.1. Pure aggregation and pure breakage Two case studies are considered with pure breakage and pure aggregation processes. It is assumed that particle growth and nucleation are negligible. A logarithmic grid is used to capture the PSD in a wide range of size domain.

ð62Þ

Note that the growth rate in Eq. (62) is size dependent. Thus with the appropriate coordinate transformation, the proposed LBM-FP method can be applied for efficient simulation of simultaneous processes with growth, nucleation, aggregation and breakage. Nucleation term is interpreted as the following boundary condition for the PBE in Eq. (60) (Ramkrishna, 2000): 1

hð0,tÞFðtÞ ¼ B0 ðtÞdðc

1

ðz0 ÞÞO ðc

ðzÞÞ:

ð63Þ

This is a Dirichlet type boundary condition and is similar to a fluid flow problem with specified inlet flow rate. This type of boundary

4.1.1. Pure breakage First, we consider a pure breakage problem. The initial distribution is taken as the following exponential distribution: nðx,0Þ ¼ expðxÞ:

ð67Þ

The breakage rate and breakage probability are GðxÞ ¼ x and bðx,x0 Þ ¼ 2=x0 , respectively, which corresponds to a binary breakage problem. The analytical solution of this problem is given by Ziff and McGrady (1985) as nðx,tÞ ¼ expðxð1 þtÞÞð1 þ tÞ2 :

ð68Þ

322

A. Majumder et al. / Chemical Engineering Science 69 (2012) 316–328

4.1.2. Pure aggregation The next problem is a pure aggregation problem with a constant aggregation kernel aðx,x0 Þ ¼ 1, where the initial PSD is taken as N0 expðx=xa Þ: xa

ð69Þ

Analytical solutions for this kernel are given by Scott (1968) as   4N0 2x0 , ð70Þ exp nðx,tÞ ¼ tþ2 xa ðt þ 2Þ2

PSD, n (#/µm)

nðx,0Þ ¼

100

Initial

10−5

LBM FP

where t ¼ N0 at and x0 ¼ x=xa . The parameters N0 ¼ 1, xa ¼1 mm are the total number of particles per unit volume and mean size of the particles, respectively. 4.1.3. Results and discussion All the problems are solved using proposed LBM-FP method and the results are compared with the pure FP (Kumar and Ramkrishna, 1996a) and FV (Qamar et al., 2009) methods. The simulations are performed using a logarithmic grid obtained by defining the transformation as O ðxÞ ¼ x, so that z ¼ logðxÞ. Final time ðt f Þ for these two problems are taken as 10 and 20 s, respectively. Final distributions for the problems obtained with N¼ 100 grid points are shown in Figs. 2 and 3. In Fig. 2, it can be seen that as the particles break, more smaller sized particles are produced. As a result, most of the particles are concentrated in the smaller size region. However, as expected, an opposite scenario is seen in Fig. 3, where the particle size increases significantly due to aggregation and the number of the smaller sized particles decreases. From these plots, it can be noted that the solutions obtained by the proposed method are in very good agreement with the analytical solutions. It is able to capture the PSD in the whole range of size domain. The L1 and L2 error norms for all the methods are summarized in Table 1. The performance of pure FP and FV methods is found to be similar in terms of accuracy for the aggregation example, while FP method is more efficient. However, the proposed method, although based on FP method for evaluation of the aggregation and breakage terms, is found to be less accurate as compared to pure FP method. This happens as LBM scheme was primarily developed for simulation of pure growth processes (i.e., advection type). For aggregation and breakage processes (without growth), the streaming step in LBM introduces undesired numerical diffusion as can be seen in the larger size region in Fig. 3. When compared with FV method, the proposed LBM-FP method is more

PSD, n (#/µm)

100 Initial LBM−FP FP

FV Analytical

10−10 10−2

100 Particle size, x (µm)

102

Fig. 3. Final distribution for pure aggregation problem at tf ¼ 20 s.

Table 1 Simulation results for pure aggregation and pure breakage processes. Cases

Parameter

LBM-FP

FP

FV

Breakage

L1-norm L2-norm CPU time (s)

3.86  10  2 9.93  10  2 0.04

3.82  10  2 8.86  10  2 0.06

7.68  10  2 0.17 0.43

Aggregation

L1-norm L2-norm CPU time (s)

1.44  10  3 3.07  10  4 1

3.18  10  4 5.58  10  5 0.06

2.57  10  4 6.25  10  5 0.37

accurate in the case of the breakage example, although it is not very obvious from Fig. 2. The reason is that LBM-FP method captures the CSD more accurately in size range where CSD is higher. This issue is further discussed in Section 4.2. In order to explore how these methods perform with varying grid size, the trade-off plots obtained using different grid sizes are shown in Figs. 4 and 5 for breakage and aggregation problems, respectively. These plots provide the useful information about the computational cost required to obtain a solution with tolerable error. In these plots, curves that are located near the bottom-left corner correspond to methods which can provide more accurate solution with less computation time. It is clear from these plots that FP method performs better than LBM-FP method for both aggregation and breakage. Thus, for pure breakage and pure aggregation no significant advantage is provided by the proposed LBM-FP method. It is shown in the next section that for growth dominant simultaneous processes e.g., processes with simultaneous growth and aggregation, the proposed method performs better than FV technique. The proposed method is also found to be advantageous for constant growth and processes with nucleation as compared to MOC-FP method (Kumar and Ramkrishna, 1997).

FV

10−5

4.2. Simultaneous growth and aggregation

Analytical

10−4

10−2 Particle size, x (µm)

100

Fig. 2. Final distribution for pure breakage problem at tf ¼ 10 s.

In this section, we consider some processes with simultaneous growth and aggregation. The various growth rates and aggregation kernels used for the case study are shown in Table 2. The initial distribution is given by Eq. (69) with N 0 ¼ 5, xa ¼ 0:01 mm for Cases 1 and 2, and N0 ¼ 1, xa ¼ 1 mm for Case 3. The analytical solutions are available for the first two cases (Ramabhadran et al., 1976). Analytical solution is not available for Case 3 and a solution obtained with refined grid is treated as a reference solution for

A. Majumder et al. / Chemical Engineering Science 69 (2012) 316–328

LBM−FP FP FV

LBM−FP FP FV

10−1 L2−norm

L1−norm

10−1

323

10−2

10−2

10−3

10−1 100 Computation time (s)

101

10−1 100 Computation time (s)

Fig. 4. Trade-off plots for pure breakage. Solid curves represent the fitted curves assuming power law relation. (a) L1-norm. (b) L2-norm.

10−3

10−4

LBM−FP FV FP

L2−norm

L1−norm

10−4 LBM−FP FV FP

10−5

10−5

100

102

10−1

Computation time (s)

100 101 Computation time (s)

Fig. 5. Trade-off plots for pure aggregation. Solid curves represent the fitted curves assuming power law relation. (a) L1-norm. (b) L2-norm.

Table 2 Parameters for simultaneous growth and aggregation processes where g ¼ 0:6 and growth exponent g ¼ 0.7. G(x)

a

L

1 2 3

1 x ð1þ gxÞg

100 x þ x0 0:1xx0

0.20 20 10

comparing the accuracy of different methods. Two different processes, i.e., growth and aggregation, are involved in these studies. These rate processes have fundamental characteristic time scales. Ramabhadran et al. (1976) defined a dimensionless variable L as the ratio of the characteristic time scales for growth and aggregation. For constant growth rate G and constant aggregation kernel a, L can be written as

L¼

G : aN0 xa

PSD, n (#/µm)

100

Case

10−5

10−10

Initial Analytical LBM−FP FV MOC−FP

10−2

100 Particle size, x (µm)

102

ð71Þ

This dimensionless number L gives us an idea of the governing rate process. For example, if L 4 1, the growth process is faster than aggregation and dominates the dynamics of the system. The values of L for the different cases are shown in Table 2. For the linear growth rate function in Case 2, the coordinate transformation is defined as dx=dz ¼ x. This transformation results in a logarithmic grid that is suited for the aggregation process. In Case 3, the transformation is defined as dx=dz ¼ ð1 þ gxÞg . With these transformations, the growth rate in the transformed coordinate becomes size-independent. On the other hand, for the constant growth rate G ¼ 1 mm s1 in Case 1, similar coordinate transformation is performed as Case 2 to obtain the logarithmic grid for which the growth rate GðzÞ ¼ expðgzÞ is size-dependent, as shown in Eq. (62).

Fig. 6. Final distribution for simultaneous growth and aggregation at tf ¼1 s, N ¼200 (Case 1).

All the case studies are solved using the proposed LBM-FP technique and the results are compared with those obtained using the FV and MOC-FP methods. Final distributions obtained with N¼ 200 grid points along with the analytical solutions are shown in Figs. 6–8. A common feature of the final distributions shown in these figures is that the distributions have a sharp front moving from the smaller sized region. This front is attributed to the growth process which results in the translation of the PSD along the computational domain with time. In addition, the presence of the aggregation process causes the shape of the PSD to stretch by increasing the number of larger particles at the cost

324

A. Majumder et al. / Chemical Engineering Science 69 (2012) 316–328

Table 3 Simulation results for simultaneous growth and aggregation. Cases

Parameter

LBM-FP

FV

MOC-FP

1

L1-norm L2-norm CPU time (s)

6.94  10  4 6.60  10  4 4.05

9.42  10  4 8.82  10  4 2.04

6.09  10  4 6.45  10  4 31.82

2

L1-norm L2-norm CPU time (s)

3.64  10  2 0.18 0.20

0.30 0.67 1.57

1.62  10  3 1.96  10  3 0.15

3

L1-norm L2-norm CPU time (s)

0.13 0.10 0.22

0.42 0.66 0.29

6.74  10  3 4.07  10  3 3.89

PSD, n (#/µm)

100

Initial Analytical

10−5

LBM−FP FV MOC−FP

10

−10

10−4

10−2 Particle size, x (µm)

100

50

PSD, n (#/µm)

40 Initial

30

Analytical LBM−FP

20

FV MOC−FP

10 0 10−4

10−3 Particle size, x (µm)

10−2

Fig. 7. Final distribution for simultaneous growth and aggregation at tf ¼ 2 s, N ¼200 (Case 2). (a) Final distribution. (b) Zoomed view.

PSD, n (#/µm)

100

10−5 Initial LBM−FP FV

10−10

10−15 −2 10

MOC−FP

100 Particle size, x (µm)

102

Fig. 8. Final distribution for simultaneous growth and aggregation at tf ¼ 1 s, N ¼200 (Case 3).

of smaller particles. It can be noted that the final distributions obtained using all the three methods match well. However, diffusive solutions are obtained for LBM-FP and FV methods in the small size region where sharp moving fronts are encountered. The diffusion is more pronounced for LBM-FP technique. The error norms and the computation times for all the case studies are presented in Table 3. We see that for Case 1, all the methods have similar accuracy while MOC-FP method is most expensive. This is due to the fact that in the MOC-FP method, the newly born particles due to aggregation, which do not fall on a grid, have to be redistributed among the neighboring grids so that some properties of interest are preserved (the first two moments of PSD in this case). Since, the grids are also moving at the characteristic velocity, this redistribution rule has to be updated frequently (Kumar and Ramkrishna, 1997). On the other hand, MOC-FP method is most accurate and efficient for Case 2 where a linear growth rate is considered. Here the redistribution rule for the FP method need not be updated since at any time t the grid points are given as xi ðtÞ ¼ xi ð0Þ expðG0 tÞ (Kumar and Ramkrishna, 1997). The LBM-FP method is found to be more accurate than FV method even though this is not apparent in Fig. 7(a) where log–log plot of the final distribution is shown. However, from Fig. 7(b) where zoomed plot is shown on a semi-log axis, it is clear that in the size range x ¼ 104 2102 mm the PSD corresponding to LBM-FP method is more close to the analytical solution than FV method. The larger deviations of the PSD due to FV method from the analytical solution in that size range dominate the error norms and thus error norms for FV method are greater than LBM-FP method. In Case 3, nonlinear growth and product kernel are considered in order to test the applicability of the proposed technique to more complex processes. Since analytical solution is not available for this example, we take the solutions obtained with N¼1000 as the reference solutions to calculate the error norm. It can be seen from Table 3 that LBM-FP is again more accurate than the FV method, while MOC-FP method is most accurate with largest computation time among the three methods. The accuracy of the MOC-FP method is because of its negligible numerical diffusion and large computation time is caused by the update of the redistribution rule for aggregation. In order to ensure that the presented results are not coincidental and to have a better idea of the performance of these methods with different number of grid points, the trade-off plots between the error norms and the computation times are shown in Figs. 9–11. From Fig. 9 it can be seen that for the aggregation dominated problem (L o 1), the trade-off curves for LBM-FP and FV schemes almost overlap each other which implies that for a given error norm, the computation time for both methods will be similar. Thus, the relatively higher computation time in LBM-FP method for the aggregation dominated process is compensated by the more accurate solution. No significant improvement is found for MOC-FP

A. Majumder et al. / Chemical Engineering Science 69 (2012) 316–328

10−3 L2−norm

L1−norm

10−3

LBM−FP FV MOC−FP

10−4 0 10

325

LBM−FP FV MOC−FP

10−4

102 Computation time (s)

101 102 Computation time (s)

103

Fig. 9. Trade-off plots for simultaneous growth and aggregation (Case 1 in Table 2). Solid curves represent the fitted curves assuming power law relation. (a) L1-norm. (b) L2-norm.

100 10−1 L2−norm

L1−norm

10−1 10−2 LBM−FP FV MOC−FP

10−3 10−1

10−2 10−3

100 101 Computation time (s)

10−1

LBM−FP FV MOC−FP

100 101 Computation time (s)

105

100

10−1

10−1 L2−norm

L1−norm

Fig. 10. Trade-off plots for simultaneous growth and aggregation (Case 2 in Table 2). Solid curves represent the fitted curves assuming power law relation. (a) L1-norm. (b) L2-norm.

10−2 10−3 10−4 −1 10

LBM−FP FV MOC−FP

100 101 102 Computation time (s)

10−2

LBM−FP FV MOC−FP

10−3

103

10−4 −1 10

100 101 102 Computation time (s)

103

Fig. 11. Trade-off plots for simultaneous growth and aggregation (Case 3 in Table 2). Solid curves represent the fitted curves assuming power law relation. (a) L1-norm. (b) L2-norm.

method with the refinement of grids for this case study. On the other hand, for the growth dominated problems (L 41) in Cases 2 and 3, it can be noted in Figs. 10 and 11 that for a given error norm, MOC-FP method is the most efficient followed by LBM-FP, while FV method is least efficient. As discussed earlier, the advantage of MOC-FP is more evident for linear growth rate. It is shown in the next section that when nucleation is present, LBM-FP performs better than the MOC-FP method as implementation of nucleation in MOC-FP method is not straight forward and requires special effort (Kumar and Ramkrishna, 1997; Lim et al., 2002). 4.3. Simultaneous growth, nucleation, aggregation and breakage In this section, it is shown that the developed method can be applied to practically relevant case of simultaneous growth,

nucleation, aggregation and breakage. We consider an example taken from Lim et al. (2002), where the initial distribution is ( 100 for 0:4 rx r0:6, nðx,0Þ ¼ ð72Þ 0:01 elsewhere: The expression for stiff nucleation as a function of time is given as nðt,0Þ ¼ 100 þ106 expð104 ðt0:215Þ2 Þ:

ð73Þ

The particles are assumed to be growing with the growth rate of G ¼ 1:0 mm s1 and the computational domain is taken as 0 rx r 2:0 mm. The parameters for aggregation and breakage are taken as aðx,x0 Þ ¼ 1:5  105 , bðx,x0 Þ ¼ 2=x0 and GðxÞ ¼ x2 . This test problem is solved using LBM-FP, FV and MOC-FP methods for final time tf ¼0.5 s with the number of grid points

326

A. Majumder et al. / Chemical Engineering Science 69 (2012) 316–328

being N ¼200. It is necessary to adaptively determine the time levels to add a grid to accommodate the newly born particles while dealing with nucleation in MOC-FP method. For implementation of nucleation, we take the approach of Lim et al. (2002) which is to add a new grid at a given time level and delete the last grid at the same time level to keep the number of grid points constant. For MOC-FP method, the absolute tolerance of 108 and relative tolerance of 104 are used in the ODE solver, since the default tolerance settings are not able to resolve the stiff nucleation. The final distributions obtained are shown in Fig. 12. The growth process causes the PSD to translate along the size axis while the aggregation and breakage processes produce larger and smaller particles, respectively. The combined effect of all these processes along with nucleation is seen in the final distribution.

Initial LBM−FP FV

105 PSD, n (#/µm)

MOC−FP

100

0

0.5

1 Particle size, x (µm)

5. Conclusions

1.5

Fig. 12. Final distribution for simultaneous growth, nucleation, aggregation and breakage after tf ¼0.5 s with N¼ 200.

Table 4 Simulation results for 1D process with simultaneous growth, nucleation, aggregation and breakage, N ¼ 200. Parameter

LBM-FP

FV

MOC-FP

L1-norm L2-norm CPU time (s)

1.13  104 4.31  104 0.68

1.98  104 7.74  104 3.99

1.44  104 7.37  104 72.28

Lattice Boltzmann method (LBM) has been introduced for solution of 1D population balance equations (PBEs) with simultaneous growth, nucleation, aggregation and breakage. Aggregation and breakage, which act as source terms in PBE, have been incorporated as force terms in the LBM formulation. These force terms can be evaluated by the techniques available in the literature for solving PBEs with aggregation and breakage. We have used fixed pivot (FP) method by Kumar and Ramkrishna (1996a) for its easy implementation and flexibility. The performance of the proposed method is compared with pure FP, finite volume (FV) and method of characteristic (MOC) combined with FP (MOC-FP) methods in terms of error norms and computation time. It has been found that FP method performs best for pure aggregation and breakage. The relatively less accurate solution for the proposed method is due to

LBM−FP FV MOC−FP

L1−norm

LBM−FP FV MOC−FP

104.8 L2−norm

105

The plateau in the final distribution seen near x ¼ 1 mm is due to the initial distribution. The spread of this square pulse is attributed to the aggregation and breakage processes. The part of the curve before this plateau is mostly owing to the newly born particles due to nucleation. However, in the final distribution the nucleation part from the FV method is slightly shifted to the left. This is likely due to the reformulation of the PBE used in FV method where the original PBE is written in mass conservation form. The analytical solution for this problem is not available. In order to compare the accuracy of the proposed technique with the FV method, we take the final distributions obtained by all the methods with a finer grid, i.e., N¼1000, as the reference solutions. The error norms and computational time for these methods are presented in Table 4 which suggest that the proposed method is the most accurate among the three methods, while it is significantly efficient. The trade-off plots showing the error norms and the computation times, which have been obtained by solving the problem with different number of grid points, are shown in Fig. 13. In this figure, the curves for the proposed LBM-FP method are situated at the bottom-left corner indicating that for a given accuracy, the computation time required for the proposed method is much lower (about an order of magnitude) than the other methods. The MOC-FP method is computationally expensive since at some time levels a new grid is added to accommodate nucleation. Consequently, the solver has to initialize the variables and evaluate the redistribution rule used in the FP method for the newly born particles due to aggregation and breakage. Thus, it can be concluded that the proposed scheme is more suitable for simulation of simultaneous processes involving growth, nucleation, aggregation and breakage due to its higher accuracy and significant advantage of computation time as compared to the FV and MOC-FP methods.

104.6

104 104.4 100 102 Computation time (s)

100

102 Computation time (s)

104

Fig. 13. Trade-off plots for simultaneous growth, nucleation, aggregation and breakage. Solid curves represent the fitted curves assuming power law relation. (a) L1-norm. (b) L2-norm.

A. Majumder et al. / Chemical Engineering Science 69 (2012) 316–328

the fact that LBM was primarily developed for growth processes (advection type) and for pure aggregation and breakage processes, the growth term is zero. For such cases, the streaming step in LBM introduces numerical diffusion. For simultaneous processes with growth and aggregation, the proposed LBM-FP method is found to have similar performance as the FV method for aggregation dominated processes, and better performance when the growth process is dominant in comparison to the aggregation process. The MOC-FP method is found to be most accurate for these simultaneous growth and aggregation processes. However, when nucleation is present, the LBM-FP method performs better than the MOC-FP method, since implementation of nucleation in this technique is difficult. The use of the MOC-FP method has some limitations as solution quality varies according to the frequency of mesh addition and computation time can be large for stiff nucleation. In addition, MOC cannot be used to solve some more complex processes, e.g., where more than one species with different growth rates are present in the system and no unique characteristic curve exists (Lim et al., 2002). Overall, the proposed LBM-FP method is found to be suitable for solution of PBEs for its wide range of applicability and efficiency. Theoretical analysis of the observation that LBM-FP has better efficiency and accuracy for growth dominant process compared to the aggregation dominant process, will be pursued in future. In addition, LBM will be extended to multi-dimensional PBEs with growth, nucleation, aggregation and breakage, where the most challenging part to be addressed is the multi-component aggregation.

Nomenclature aðx,x0 Þ rate of aggregation of particles with size x and x0 , s  1 bðx,x0 Þ probability of formation of particles with size x through breakage of particles with size x0 , s  1 B0 rate of nucleation, # s  1 c nonzero velocity component of the particles along the coordinate axes, mm s  1 cs speed of sound, mm s  1 D diffusivity, mm2 s1 fi discrete Boltzmann distribution function, # mm1 eq equilibrium discrete Boltzmann distribution, # mm1 fi ^f re-defined discrete distribution function, # mm1 i

G h Ma n N P Q r t u v x z

growth rate of particles, mm s1 distribution function in transformed coordinate system, # mm1 Mach number, dimensionless particle size distribution, # mm1 number of grid points, dimensionless pressure, kg m  1 s  2 various rate processes, e.g., aggregation, #mm1 s1 space coordinate, mm time, s local average velocity, mm s1 advection velocity, analogous to G in PBE, mm s1 characteristic size of particles, mm size in the transformed coordinate system, mm

Greek letters

b L F

r

relaxation parameter, dimensionless ratio of the characteristic time scales for growth and aggregation, dimensionless time-dependent part of growth rate, mm s1 concentration of the transported species, g mm3

t O

327

relaxation time, s size-dependent part of growth rate, mm s1

References Aamir, E., Nagy, Z.K., Rielly, C.D., Kleinert, T., Judat, B., 2009. Combined quadrature method of moments and method of characteristics approach for efficient solution of population balance models for dynamic modeling and crystal size distribution control of crystallization processes. Ind. Eng. Chem. Res. 48 (18), 8575–8584. Abegg, C.F., Stevens, J.D., Larson, M.A., 1968. Crystal size distributions in continuous crystallizers when growth rate is size dependent. AIChE J. 14 (1), 118–122. Aidun, C.K., Clausen, J.R., 2009. Lattice-Boltzmann method for complex flows. Annu. Rev. Fluid Mech. 42 (1), 439–472. Alexopoulos, A.H., Roussos, A., Kiparissides, C., 2009. Part V: dynamic evolution of the multivariate particle size distribution undergoing combined particle growth and aggregation. Chem. Eng. Sci. 64 (14), 3260–3269. Ansumali, S., Karlin, I.V., 2000. Stabilization of the lattice Boltzmann method by the h theorem: a numerical test. Phys. Rev. E 62 (6), 7999. Ansumali, S., Karlin, I.V., 2002a. Entropy function approach to the lattice Boltzmann method. J. Stat. Phys. 107, 291–308. Ansumali, S., Karlin, I.V., 2002b. Single relaxation time model for entropic lattice Boltzmann methods. Phys. Rev. E 65, 056312. Ansumali, S., Karlin, I.V., Arcidiacono, S., Abbas, A., Prasianakis, N.I., 2007. Hydrodynamics beyond Navier–Stokes: exact solution to the lattice Boltzmann hierarchy. Phys. Rev. Lett. 98 (12), 124502. Batterham, R.J., Hall, J.S., Barton, G., 1981. Pelletizing kinetics and simulation of full scale balling circuits. In: Proceedings of the 3rd International Symposium on Agglomeration, Nurnberg, W. Germany, p. A136. Benzi, R., Succi, S., Vergassola, M., 1992. The lattice Boltzmann equation: theory and applications. Phys. Rep. 222 (3), 145–197. Bhatnagar, P.L., Gross, E.P., Krook, M., 1954. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94 (3), 511. Bilgili, E., Scarlett, B., 2005. Population balance modeling of non-linear effects in milling processes. Powder Technol. 153 (1), 59–71. Chen, H., Kandasamy, S., Orszag, S., Shock, R., Succi, S., Yakhot, V., 2003. Extended Boltzmann kinetic equation for turbulent flows. Science 301 (5633), 633–636. Chen, S., Doolen, G.D., 1998. Lattice Boltzmann method for fluid flows. Annu. Rev. Fluid Mech. 30 (1), 329–364. Chikatamarla, S.S., Babu, V., 2005. Entropic lattice Boltzmann method on nonuniform grids. Computational Science—ICCS 2005. Lecture Notes in Computer Science, vol. 3516. , Springer, Berlin/Heidelberg, pp. 72–79. Costa, C.B.B., Maciel, M.R.W., Filho, R.M., 2007. Considerations on the crystallization modeling: population balance solution. Comput. Chem. Eng. 31 (3), 206–218. Dawson, S.P., Chen, S., Doolen, G.D., 1993. Lattice Boltzmann computations for reaction-diffusion equations. J. Chem. Phys. 98 (2), 1514–1523. Filbet, F., Laurencot, P., 2004. Numerical simulation of the Smoluchowski coagulation equation. SIAM J. Sci. Comput. 25 (6), 2004–2028. Frisch, U., Hasslacher, B., Pomeau, Y., 1986. Lattice-gas automata for the Navier– Stokes equation. Phys. Rev. Lett. 56 (14), 1505. Gimbun, J., Nagy, Z.K., Rielly, C.D., 2009. Simultaneous quadrature method of moments for the solution of population balance equations, using a differential algebraic equation framework. Ind. Eng. Chem. Res. 48 (16), 7798–7812. Gunawan, R., Fusman, I., Braatz, R.D., 2004. High resolution algorithms for multidimensional population balance equations. AIChE J. 50 (11), 2738–2749. He, X., Luo, L.-S., Dembo, M., 1996. Some progress in lattice Boltzmann method. Part I. Nonuniform mesh grids. J. Comput. Phys. 129 (2), 357–363. Higuera, F., Succi, S., Benzi, R., 1989. Lattice gas-dynamics with enhanced collisions. Europhys. Lett. 9, 345–349. Immanuel, C.D., Cordeiro, C.F., Sundaram, S.S., Meadows, E.S., Crowley, T.J., Doyle III, F.J., 2002. Modeling of particle size distribution in emulsion co-polymerization: comparison with experimental data and parametric sensitivity studies. Comput. Chem. Eng. 26 (7–8), 1133–1152. Immanuel, C.D., Doyle III, F.J., 2003. Computationally efficient solution of population balance models incorporating nucleation, growth and coagulation: application to emulsion polymerization. Chem. Eng. Sci. 58 (16), 3681–3698. Immanuel, C.D., Doyle III, F.J., 2005. Solution technique for a multi-dimensional population balance model describing granulation processes. Powder Technol. 156 (2–3), 213–225. Kariwala, V., Cao, Y., Nagy, Z.K. Automatic differentiation based quadrature method of moments for solving population balance equations. AIChE J., doi:10.1002/aic.12613, In press. Karlin, I.V., Ansumali, S., Frouzakis, C.E., Chikatamarla, S.S., 2006. Elements of the lattice Boltzmann method. I. Linear advection equation. Commun. Comput. Phys. 1 (4), 616–655. ¨ ttinger, H.C., 1999. Perfect entropy functions of the Karlin, I.V., Ferrante, A., O lattice Boltzmann method. Europhys. Lett. 47 (2), 182–188. Kim, S.H., Pitsch, H., Boyd, I.D., 2008. Accuracy of higher-order lattice Boltzmann methods for microscale flows with finite Knudsen numbers. J. Comput. Phys. 227 (19), 8655–8671.

328

A. Majumder et al. / Chemical Engineering Science 69 (2012) 316–328

Kumar, J., Peglow, M., Warnecke, G., Heinrich, S., 2008a. The cell average technique for solving multi-dimensional aggregation population balance equations. Comput. Chem. Eng. 32 (8), 1810–1830. Kumar, J., Peglow, M., Warnecke, G., Heinrich, S., 2008b. An efficient numerical technique for solving population balance equation involving aggregation, breakage, growth and nucleation. Powder Technol. 182 (1), 81–104. Kumar, J., Peglow, M., Warnecke, G., Heinrich, S., Morl, L., 2006. Improved accuracy and convergence of discretized population balance for aggregation: the cell average technique. Chem. Eng. Sci. 61 (10), 3327–3342. Kumar, J., Warnecke, G., Peglow, M., Heinrich, S., 2009. Comparison of numerical methods for solving population balance equations incorporating aggregation and breakage. Powder Technol. 189 (2), 218–229. Kumar, S., Ramkrishna, D., 1996a. On the solution of population balance equations by discretization. I. A fixed pivot technique. Chem. Eng. Sci. 51 (8), 1311–1332. Kumar, S., Ramkrishna, D., 1996b. On the solution of population balance equations by discretization. II. A moving pivot technique. Chem. Eng. Sci. 51 (8), 1333–1342. 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 (24), 4659–4679. Lang, Y.D., Cervantes, A.M., Biegler, L.T., 1999. Dynamic optimization of a batch cooling crystallization process. Ind. Eng. Chem. Res. 38 (4), 1469–1477. LeVeque, R.J., 2002. Finite-Volume Methods for Hyperbolic Problems. Cambridge University Press, Cambridge, UK. Lim, Y., Le Lann, J.-M., Meyer, X.M., Joulia, X., Lee, G., Yoon, E.S., 2002. On the solution of population balance equations (PBE) with accurate front tracking methods in practical crystallization processes. Chem. Eng. Sci. 57 (17), 3715–3732. Mahoney, A.W., Ramkrishna, D., 2002. Efficient solution of population balance equations with discontinuities by finite elements. Chem. Eng. Sci. 57 (7), 1107–1119. Majumder, A., Kariwala, V., Ansumali, S., Rajendran, A., 2010a. Entropic lattice Boltzmann method for crystallization processes. Chem. Eng. Sci. 65 (13), 3928–3936. Majumder, A., Kariwala, V., Ansumali, S., Rajendran, A., 2010b. Fast high-resolution method for solving multidimensional population balances in crystallization. Ind. Eng. Chem. Res. 49 (8), 3862–3872. Majumder, A., Kariwala, V., Ansumali, S., Rajendran, A. Lattice Boltzmann method for multi-dimensional population balance models in crystallization. Chem. Eng. Sci., doi:10.1016/j.ces.2011.04.041, In press. Mantle, M.D., Sederman, A.J., Gladden, L.F., 2001. Single- and two-phase flow in fixed-bed reactors: MRI flow visualisation and lattice-Boltzmann simulations. Chem. Eng. Sci. 56 (2), 523–529. Mantzaris, N.V., Daoutidis, P., Srienc, F., 2001. Numerical solution of multi-variable cell population balance models. I. Finite difference methods. Comput. Chem. Eng. 25 (11–12), 1411–1440. Marchisio, D.L., Fox, R.O., 2005. Solution of population balance using the direct quadrature method of moments. J. Aerosol Sci. 36 (1), 43–73. Marchisio, D.L., Vigil, R.D., Fox, R.O., 2003. Quadrature method of moments for aggregation-breakage processes. J. Colloid Interface Sci. 258 (2), 322–334. McNamara, G.R., Zanetti, G., 1988. Use of the Boltzmann equation to simulate lattice-gas automata. Phys. Rev. Lett. 61 (20), 2332. Meimaroglou, D., Kiparissides, C., 2007. Monte Carlo simulation for the solution of the bi-variate dynamic population balance equation in batch particulate systems. Chem. Eng. Sci. 62 (18–20), 5295–5299. Motz, S., Mitrovic´, A., Gilles, E.D., 2002. Comparison of numerical methods for the simulation of dispersed phase systems. Chem. Eng. Sci. 57 (20), 4329–4344. Myerson, A.S., 2002. Handbook of Industrial Crystallization, second ed. ButterworthHeinemann, Boston, USA.

Nagy, Z.K., 2009. Model based robust control approach for batch crystallization product design. Comput. Chem. Eng. 33 (10), 1685–1691. Nicmanis, M., Hounslow, M.J., 1998. Finite-element methods for steady-state population balance equations. AIChE J. 44 (10), 2258–2272. Nopens, I., Beheydt, D., Vanrolleghem, P.A., 2005. Comparison and pitfalls of different discretised solution methods for population balance models: a simulation study. Comput. Chem. Eng. 29 (2), 367–377. Pinto, M.A., Immanuel, C.D., Doyle III, F.J., 2007. A feasible solution technique for higher-dimensional population balance models. Comput. Chem. Eng. 31 (10), 1242–1256. Pinto, M.A., Immanuel, C.D., Doyle III, F.J., 2008. A two-level discretisation algorithm for the efficient solution of higher-dimensional population balance models. Chem. Eng. Sci. 63 (5), 1304–1314. Poon, J.M.H., Ramachandran, R., Sanders, C.F.W., Glaser, T., Immanuel, C.D., Doyle III, F.J., Litster, J.D., Stepanek, F., Wang, F.-Y., Cameron, I.T., 2009. Experimental validation studies on a multi-dimensional and multi-scale population balance model of batch granulation. Chem. Eng. Sci. 64 (4), 775–786. Qamar, S., Mukhtar, S., Ali, Q., Seidel-Morgenstern, A., 2011. A Gaussian quadrature method for solving batch crystallization models. AIChE J. 57 (1), 149–159. Qamar, S., Warnecke, G., 2007. Numerical solution of population balance equations for nucleation, growth and aggregation processes. Comput. Chem. Eng. 31 (12), 1576–1589. Qamar, S., Warnecke, G., Elsner, M.P., 2009. On the solution of population balances for nucleation, growth, aggregation and breakage processes. Chem. Eng. Sci. 64 (9), 2088–2095. Ramabhadran, T.E., Peterson, T.W., Seinfeld, J.H., 1976. Dynamics of aerosol coagulation and condensation. AIChE J. 22 (5), 840–851. Ramkrishna, D., 2000. Population Balances: Theory and Applications to Particulate Systems in Engineering. Academic Press, San Diego, CA. Randolph, A.D., Larson, M.A., 1962. Transient and steady state size distributions in continuous mixed suspension crystallizers. AIChE J. 8 (5), 639–645. Roussas, G., 2003. An Introduction to Probability and Statistical Inference. Academic Press, San Diego, CA. Scott, W.T., 1968. Analytic studies of cloud droplet coalescence I. J. Atmos. Sci. 25 (1), 54–65. Silva, L.F.L.R., Rodrigues, R.C., Mitre, J.F., Lage, P.L.C., 2010. Comparison of the accuracy and performance of quadrature-based methods for population balance problems with simultaneous breakage and aggregation. Comput. Chem. Eng. 34 (3), 286–297. Singh, S., Subramanian, G., Ansumali, S., 2011. A lattice Boltzmann method for dilute polymer solutions. Philos. Trans. R. Soc. A 369, 2301–2310. Smith, M., Matsoukas, T., 1998. Constant-number Monte Carlo simulation of population balances. Chem. Eng. Sci. 53 (9), 1777–1786. Succi, S., 2001. Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Oxford University Press, New York. van Peborgh Gooch, J.R., Hounslow, M.J., 1996. Monte Carlo simulation of sizeenlargement mechanisms in crystallization. AIChE J. 42 (7), 1864–1874. Woo, X.Y., Tan, R.B.H., Braatz, R.D., 2009. Modeling and computational fluid dynamics – population balance equation – micromixing simulation of impinging jet crystallizers. Cryst. Growth Des. 9 (1), 156–164. Yudistiawan, W.P., Kwak, S.K., Patil, D.V., Ansumali, S., 2010. Higher-order Galilean-invariant lattice Boltzmann model for microflows: single-component gas. Phys. Rev. E 82 (4), 046701. Zhao, H., Maisels, A., Matsoukas, T., Zheng, C., 2007. Analysis of four Monte Carlo methods for the solution of population balances in dispersed systems. Powder Technol. 173 (1), 38–50. Ziff, R.M., McGrady, E.D., 1985. The kinetics of cluster fragmentation and depolymerisation. J. Phys. A: Math. Gen. 18 (15), 3027.