Solution of Bivariate Population Balance Equations Using the Finite

262

Ind. Eng. Chem. Res. 2009, 48, 262–273

Solution of Bivariate Population Balance Equations Using the Finite Size Domain Complete Set of Trial Functions Method of Moments (FCMOM) Matteo Strumendo* and Hamid Arastoopour Department of Chemical and Biological Engineering, Illinois Institute of Technology, Chicago, Illinois 60616

The FCMOM (Finite size domain Complete set of trial functions Method Of Moments) is an efficient and accurate numerical technique to solve PBE (population balance equations) and was validated for monovariate PBE [Strumendo, M.; Arastoopour, H. Solution of PBE by MOM in Finite Size Domains. Chem. Eng. Sci. 2008, 63 (10), 2624]. In the present paper, the FCMOM is extended and used to solve bivariate PBE. In the FCMOM, the method of moments is formulated in a finite domain of the internal coordinates and the particle distribution function is represented as a truncated series expansion by a complete system of orthonormal functions. In the extension to bivariate PBE, the capabilities of the FCMOM are maintained, particularly as far as the algorithm efficiency and the accuracy in the bivariate particle distribution function reconstruction. The FCMOM was validated with the following bivariate applications: particle growth, particle dissolution, particle aggregation, and simultaneous aggregation and growth. 1. Introduction Population balance equations (PBE) are used to model and to simulate a very wide range of applications in chemical engineering, from fluidized beds to droplet dispersions, from cell population evolution to crystallization, from size reduction to powder aggregation processes, and from polymerization reactions to aerosol formation.1 Particularly, bivariate PBE can be formulated and solved to investigate applications such as the sintering and aggregation of nanoparticles, the aggregation of a bicomponent aerosol, phenomena of anisotropic crystal growth, crystallization processes in which both the particle size and shape must be monitored, and fluidization of a population of particles varying in size and density. The bivariate formulation of the PBE is required when one internal coordinate (typically, the size) is not sufficient to define the state of a particle in the population, and the description of the state of a particle is given in terms of two coordinates (i.e., particle size and surface area, or two characteristic particle lengths, or particle size and density). A comprehensive review of the subject of PBE (formulation of PBE, solution techniques, applications, theoretical considerations) is contained in a book by Ramkrishna.1 This paper, however, is focused on the solution techniques of bivariate PBE and, more specifically, on those numerical techniques based on the method of moments (MOM) which demand a small computational effort. Such a requirement is important from the perspective of coupling the PBE with the transport equations (mass, momentum, energy, and species) of the fluid and particulate phases. In recent years, different variations of the MOM were proposed to solve PBE and a detailed discussion on this subject (mainly referring to the monovariate case) is contained in an article by Strumendo and Arastoopour.2 Focusing on the applications of the MOM to bivariate PBE, two variations of the MOM, the quadrature method of moments (QMOM) and the direct quadrature method of moments (DQMOM), were used to solve bivariate (and multivariate) PBE. * To whom correspondence should be addressed. Tel.: (312) 5673066. Fax: (312) 567-8874. E-mail: [email protected].

The QMOM was proposed by McGraw3 and then was validated (for monovariate PBE) by Barrett and Webb,4 by Marchisio et al.,5 and by Marchisio et al.6 In the QMOM, no explicit assumption is made regarding the shape of the size distribution function and the moments evolution is computed correctly and efficiently in many practical applications. The QMOM was extended for solving bivariate (and multivariate) PBE by Wright et al.,7 McGraw and Wright,8 and Yoon and McGraw9,10 using different approaches. McGraw and Wright8 developed an approach based on the JMT (Jacobian matrix transformation), which, however, is limited to internally mixed aerosols. A bivariate population balance equation was formulated by Wright et al.7 to describe simultaneous particle coagulation and sintering using the particle volume and the particle surface area as internal variables. The procedure used by Wright et al.7 to determine the abscissas and the weights in the quadratures is more complex than in the monovariate case and requires the solution of a nonlinear system at each time step; the numerical solution obtained by the bivariate formulation of the QMOM was shown to be consistent with “exact” solutions obtained by a discrete model and by Monte Carlo simulations. Yoon and McGraw9,10 proposed an alternative approach to extend the QMOM for solving bivariate and multivariate PBE using the principal components analysis (PCA). In the PCAQMOM, the quadrature points are optimally assigned through the PCA and back projection. Thus, the assignment of the quadrature abscissas and weights in the PCA-QMOM extension results computationally easier than in the extension proposed by Wright et al.7 Two disadvantages of the QMOM are related to (1) the particle distribution reconstruction and (2) the complexity and efficiency of the algorithm (in bivariate and multivariate applications). In fact, concerning the particle distribution reconstruction, in the QMOM, the solution is given in terms of the moments, while the particle distribution function disappears from the governing equations. The reconstruction of the distribution function from the moments11-13 can be a complex problem and is not always possible. Additionally, as pointed out by Marchisio and Fox,14 the QMOM algorithm,

10.1021/ie800272a CCC: $40.75  2009 American Chemical Society Published on Web 11/08/2008

Ind. Eng. Chem. Res., Vol. 48, No. 1, 2009 263

which is simple and efficient in monovariate applications, becomes quite complex in bivariate and multivariate applications requiring the solution of a nonlinear system or the use of the PCA-QMOM. The DQMOM, proposed by Marchisio and Fox,14 differs from the QMOM because the abscissas and the weights of the quadrature approximation are tracked directly (rather than the moments, as in the QMOM) and because an explicit expression for the particle distribution function is given in terms of a summation of Dirac delta functions. In this way, the DQMOM algorithm maintains its features of simplicity and efficiency when applied to bivariate and multivariate PBE;14 moreover, through the DQMOM, it is possible to take into account the dependency of the particle velocities on the internal variables. However, again, the DQMOM is an efficient method to compute the moments (tracking the abscissas and weights of the nodes) but it does not provide directly an accurate particle distribution reconstruction. Additionally, Marchisio and Fox14 pointed out that the DQMOM can yield unrealizable abscissas in bivariate and multivariate applications and suggested that the question of realizability (and its connection with the choice of the moments) must be carefully considered for each new application. Diemer and Olson13 investigated the problem of the reconstruction of the bivariate particle distribution from the moments for simultaneous coagulation, coalescence, and breakup. In the reconstruction of the bivariate particle distribution function, Diemer and Olson13 first reconstructed the marginal distribution in the volume variable, using the results obtained for the monovariate case by Diemer and Olson,11,12 and then assumed that the conditional area distribution can be expressed in terms of the beta distribution. Diemer and Olson13 were able to show reconstructed steadystate bivariate distributions (for size-independent kernels); however, the procedure followed in the bivariate particle distribution reconstruction is rather complex and relies on several restrictive assumptions. In this paper, the FCMOM (Finite size domain Complete set of trial functions Method Of Moments) is extended for solving bivariate PBE. The FCMOM was proposed and validated for monovariate PBE by Strumendo and Arastoopour.2 The solution of the PBE obtained through the FCMOM provides both the moments and the reconstructed particle size distribution. The particle size distribution is reconstructed from the moments through a simple relation, in both transient and steady-state conditions. In the extension to bivariate PBE, the FCMOM maintains its features of accuracy in the reconstruction of the bivariate particle distribution and of computational efficiency. In addition, in the FCMOM, the internal variables’ domains are defined by bounded (finite) regions, thus avoiding having the internal variables assume unrealistic values. In section 2, the general features of the bivariate formulation of the FCMOM are discussed. In sections 3, 4, and 5, some bivariate applications are developed to validate the FCMOM. In section 3, the FCMOM is applied to (1) particle growth (constant, linear, diffusion-controlled growth rate) and (2) particle dissolution. In section 4, the method is applied to particle aggregation. In section 5, an application to simultaneous particle aggregation and growth is discussed. In the present paper, spatial uniformity is assumed (homogeneous processes); therefore, the only external coordinate is the time t.

Additionally, although the formulation of the governing equations presented in section 2 is general, the domains of the internal variables considered in the applications are relatively simple (rectangular domains, in which the internal variables at the boundaries do not depend on each other); an exception is discussed in section 3.3 regarding the diffusion-controlled growth. Work is in progress to apply the FCMOM also to more complex domains. 2. Extension of the FCMOM for Solving Bivariate PBE In this section, the FCMOM, proposed and validated for monovariate PBE by Strumendo and Arastoopour,2 is extended to bivariate PBE. ξ1 and ξ2 are the internal variables,15 t is the time, and f(ξ1,ξ2,t) is the bivariate particle distribution. ξ1,min, ξ2,min and ξ1,max, ξ2,max are the minimum and maximum values of the internal variables, respectively. The internal variables ξ1 and ξ2 can be expressed in dimensionless form using the following coordinate transformations: ξ1,max(ξ2, t) + ξ1,min(ξ2, t) 2 ξ¯ 1 ) ξ1,max(ξ2, t) - ξ1,min(ξ2, t) 2 ξ2,max(ξ1, t) + ξ2,min(ξ1, t) ξ2 2 ξ¯2 ) ξ2,max(ξ1, t) - ξ2,min(ξ1, t) 2 ξ1 -

(1)

Each of the dimensionless internal variables ξ1 and ξ2 can assume values between -1 and 1. The bivariate distribution function can be defined in dimensionless terms as f′(ξ1,ξ2,t) ) f(ξ1,ξ2,t)/fsc, where fsc is an appropriate scale factor. In the bivariate formulation of the FCMOM, the dimensionless bivariate distribution function is expressed as M-1 m

f′(ξ1, ξ2, t) =

∑ ∑ φ (ξ ) φ n

1

m-n(ξ2)

cmn(t)

(2)

m)0 n)0

In eq 2, φ are the orthonormal functions associated with the Legendre polynomials2 and cmn are the coefficients of the series (2). The bivariate dimensionless moments are defined as µi1,i2 )

∫ ∫ 1

1

-1

-1

f′(ξ1)i1(ξ2)i2 dξ1 dξ2

(3)

The coefficients cmn can be expressed in terms of the moments by cmn )

2n 2+ 1 21 2(m -2n) + 1 2 1 n

(m-n)

n

∑ (-1)

n-ν

ν)0 (m-n)

∑

σ)0

×

1 (2ν)! × [(2ν - n)!] [(n - ν)!][(ν)!]

(-1)m-n-σ

(2σ)! × {[2σ - (m - n)]!}

1 µ (4) {[(m - n) - σ]!}[(σ)!] 2ν-n,2σ-(m-n) In eq 4, terms with moments of negative order are to be omitted. As in the monovariate case,2 the bivariate particle distribution f ′ is completely determined by the coefficients cmn or, alternatively, by the moments µi1,i2.

264 Ind. Eng. Chem. Res., Vol. 48, No. 1, 2009

A general formulation of the bivariate PBE for spatially homogeneous processes is15 ∂f(ξ1, ξ2, t) ∂G1(ξ1, ξ2, t) f(ξ1, ξ2, t) ∂G2(ξ1, ξ2, t) f(ξ1, ξ2, t) + + ) B(ξ1, ξ2, t) - D(ξ1, ξ2, t) ∂t ∂ξ1 ∂ξ2

(5)

where G1 and G2 are the particle growth rates; B and D are the birth and death terms due to nucleation, breakage, aggregation, etc. The moments evolution equations are derived from eq 5: dµi1,i2 dt

∫ ∫

1 fsc

+ MB1 + MB2 + IG1 + IG2 + MB12 + MB21 )

1

1

-1

-1

(B′ - D′)(ξ1)i1(ξ2)i2 dξ1 dξ2

(6)

where MB1 ) -

(

)

f′(1, ξ2, t) - f′(-1, ξ2, t)(-1)i1 ∂ξ1,max + ξ1,min (ξ2)i2 dξ2 + i1 -1 ξ1,max - ξ1,min ∂t

∫

1

(

)

∂ξ1,max + ξ1,min dξ2 ∂t

∫

1

f′(ξ1)i1-1 dξ1 -1

∫

(

1

(

)

∂ξ2,max + ξ2,min dξ1 ∂t

∫

1

f′(ξ2)i2-1 dξ2 -1

)

∫ f′(ξ )

) )

∫

(ξ2)i2 ∂ξ1,max - ξ1,min dξ2 -1 ξ ξ ∂t 1,max 1,min

∫

1

(ξ1)i1 × -1 (ξ 2,max - ξ2,min)

∫

( (

1

(ξ1)i1 ∂ξ2,max - ξ2,min dξ1 -1 ξ ∂t 2,max - ξ2,min

∫

1

1

∫

1

-1

2(ξ1)i1 [G2′(ξ1, 1, t) f′(ξ1, 1, t) - (-1)i2G2′(ξ1, -1, t) f′(ξ1, -1, t)] dξ1 - i2 -1 ξ 2,max - ξ2,min

∫

1

∫

1

-1

dξ1 (6a)

1

-1

f′(ξ2)i2 dξ2 (6b)

∫

1

-1

G1′f′(ξ1)i1-1 dξ1 (6c)

2 × ξ2,max - ξ2,min (ξ1)i1 dξ1

MB12 )

i1

1

2 × ξ1,max - ξ1,min (ξ2)i2 dξ2

IG2 )

1

-1

1

2(ξ2)i2 [G1′(1, ξ2, t) f′(1, ξ2, t) - (-1)i1 G1′(-1, ξ2, t) f′(-1, ξ2, t)] dξ2 - i1 -1 ξ 1,max - ξ1,min

∫

) )

f′(ξ1, 1, t) - f′(ξ1, -1, t)(-1)i2+1 ∂ξ2,max - ξ2,min (ξ1)i1 dξ1 + -1 ξ2,max - ξ2,min ∂t

∫

(i2 + 1) IG1 )

( (

f′(1, ξ2, t) - f′(-1, ξ2, t)(-1)i1+1 ∂ξ1,max - ξ1,min (ξ2)i2 dξ2 + -1 ξ1,max - ξ1,min ∂t

f′(ξ1, 1, t) - f′(ξ1, -1, t)(-1)i2 ∂ξ2,max + ξ2,min (ξ1)i1 dξ1 + i2 -1 ξ2,max - ξ2,min ∂t

∫

1

1

(i1 + 1) MB2 ) -

(ξ2)i2 × -1 ξ 1,max - ξ1,min

∫

∫

1

-1

G2′f′(ξ2)i2-1 dξ2 (6d)

G2′(-1, ξ2, t) f ′(-1, ξ2, t)(-1)i1+1 - G2′(1, ξ2, t) f ′(1, ξ2, t) ∂ξ1,max - ξ1,min × -1 ξ1,max - ξ1,min ∂ξ2 ∂ξ1,max - ξ1,min 1 1 1 (ξ2)i2 dξ2 -1 G2′f′(ξ1)i1 dξ1 + (ξ2)i2 dξ2 + (i1 + 1) -1 ξ1,max - ξ1,min ∂ξ2

∫

1

∫

∫

G2′(-1, ξ2, t) f ′(-1, ξ2, t)(-1)i1 - G2′(1, ξ2, t) f′(1, ξ2, t) ∂ξ1,max + ξ1,min × -1 ξ1,max - ξ1,min ∂ξ2 ∂ξ1,max + ξ1,min 1 1 (ξ2)i2 dξ2 (ξ2)i2 dξ2 + i1 -1 ξ1,max - ξ1,min ∂ξ2

∫

1

∫

∫

1

-1

G2′f′(ξ1)i1-1 dξ1 (6e)

G1′(ξ1, -1, t) f ′(ξ1, -1, t)(-1)i2+1 - G1′(ξ1, 1, t) f ′(ξ1, 1, t) ∂ξ2,max - ξ2,min MB21 ) -1 × ξ2,max - ξ2,min ∂ξ1 ∂ξ2,max - ξ2,min 1 1 1 (ξ1)i1 dξ1 -1 G1′f′(ξ2)i2 dξ2 + (ξ1)i1 dξ1 + (i2 + 1) -1 ξ2,max - ξ2,min ∂ξ1

∫

1

∫

∫

G1′(ξ1, -1, t) f′(ξ1, -1, t)(-1)i2 - G1′(ξ1, 1, t) f′(ξ1, 1, t) ∂ξ2,max + ξ2,min × -1 ξ2,max - ξ2,min ∂ξ1 ∂ξ2,max + ξ2,min 1 1 (ξ1)i1 dξ1 (ξ1)i1 dξ1 + i2 -1 ξ2,max - ξ2,min ∂ξ1

∫

1

∫

∫

1

-1

G1′f′(ξ2)i2-1 dξ2 (6f)

In eq 6, the quantities G1′, G2′, B′, and D′ are defined such that G1′(ξ1,ξ2,t) ) G1(ξ1,ξ2,t), G2′(ξ1,ξ2,t) ) G2(ξ1,ξ2,t), B′(ξ1,ξ2,t) ) B(ξ1,ξ2,t), and D′(ξ1,ξ2,t) ) D(ξ1,ξ2,t).


Figure 1. Bivariate distribution of particle characteristic lengths vs particle characteristic lengths: initial condition (circles) and numerical solution at the final time by the FCMOM (solid line) when the growth rate is diffusion-controlled along L1 and constant along L2. The final time is 500 s. M ) 12.

Figure 2. Bivariate distribution of particle characteristic lengths vs particle characteristic lengths: numerical solution at the final time by the FCMOM (solid line) and exact solution (circles) when the growth rate is diffusion-controlled along L1 and constant along L2. The final time is 500 s. M ) 12.

The structure of eq 6 is as follows. On the left-hand side, there are (1) the moments accumulation rate; (2) the terms due to the coordinate transformation (moving boundaries) MB1 and MB2, along the first and second coordinates, respectively; (3) the contributions from the integration of the particle growth rate terms IG1 and IG2, along the first and second coordinates, respectively; and (4) the terms MB12 and MB21 due to the mutual dependence of the internal variables’ boundaries (these last terms do not appear in the monovariate formulation of the FCMOM). On the right-hand side, there are the terms related to the birth or death of particles caused by nucleation, aggregation, breakage, etc.

Moments evolution equations must be coupled with the moving boundary conditions for ξ1,min, ξ2,min, ξ1,max, and ξ2,max. These are presented in sections 3, 4, and 5 when discussing each specific application. The internal variables’ domains considered in the applications of sections 3, 4, and 5 (except the application to the diffusioncontrolled growth in section 3.3) are rectangular domains, in which the internal variables at the boundaries do not depend on each other. For such domains, the terms MB12 and MB21 in eq 6 are neglected and the moments evolution equation (6) coupled with the moving boundary conditions form a set of


ordinary differential equations, whose variables are the moments µi1,i2 and the boundaries ξ1,min, ξ2,min, ξ1,max, and ξ2,max. The initial conditions for the moments are computed from the initial bivariate particle distribution function by eq 3. 3. Applications to Particle Growth and Particle Dissolution In this section, the FCMOM is applied to specific cases of pure particle growth (subsections 3.1, 3.2, and 3.3) and of pure particle dissolution (subsection 3.4). In some applications, the particles considered are defined by two geometric quantities, typically two characteristic lengths (subsections 3.1, 3.3, and 3.4). In other applications each particle is composed of two chemical species; therefore, the volume (or the mass) of the two components defines the state of each particle (subsection 3.2). Such quantities (the characteristic lengths or the volume of each of the two components) can grow or decrease (negative growth) with different rates (anisotropic growth or dissolution). For monovariate PBE, it is well-known that, when the growth rate is constant or linear, the moments equations are closed, while the diffusion-controlled growth is a relevant example in which the moments equations are not closed.2,3 Similar results are obtained in the bivariate case; subsections 3.1 and 3.2 consider the constant and linear growth rate and in subsection 3.3 an application involving the diffusion-controlled growth is discussed. In subsections 3.1, 3.3, and 3.4, the internal variables ξ1 and ξ2 are identified with two particle characteristic lengths L1 and L2; therefore, L1 ) ξ1 and L2 ) ξ2. Accordingly, the particle distribution f(ξ1,ξ2,t) is identified with the distribution of the particle characteristic lengths fL(L1,L2,t). In subsection 3.2, the internal variables ξ1 and ξ2 are identified with the volumes of each of the two components in the particle V1 and V2; therefore, V1 ) ξ1 and V2 ) ξ2. Accordingly, the particle distribution f(ξ1,ξ2,t) is identified with the particle composition distribution fV(V1,V2,t). 3.1. Constant Growth. Puel et al.16 emphasized that in crystallization a bivariate description is required when dealing with anisotropic crystals (i.e., organic products). For example, the shape factor in rodlike particles nucleating and growing in a crystallizer can be predicted by a bivariate formulation of the PBE, in which appropriate nucleation and growth rates are used; in the bivariate PBE, different growth laws hold along different particle characteristic lengths.16,17 Particularly, a constant growth rate can be assumed in crystallization when the growth is integration-controlled. Constant growth rate laws are expressed as G1 )

dL1 ) K1 dt

dL2 ) K2 dt The following moments equation results from eq 6: G2 )

dµi1,i2 dt

)0

(7) (8)

(9)

Therefore, in the dimensionless L1-L2 space, the moments do not change vs time and the bivariate particle distribution fL′(L1,L2,t) is also constant. The moving boundary conditions are given by eqs 7 and 8 applied to L1,min, L1,max and to L2,min, L2,max, respectively. Both L1,min, L1,max and L2,min, L2,max move with a constant, although different, speed. Therefore, the initial bivariate particle distribution function is convected in the L1-L2 space, toward increasing size values, at

a constant speed whose components along L1 and L2 are given by K1 and K2, respectively. 3.2. Linear Growth. Linear growth rate laws occur in systems in which the particle growth is due to a volume chemical reaction.18 In a bicomponent (or multicomponent) aerosol, each component volume (or mass) changes following a different growth law.19 Linear growth rate laws for bicomponent particles are expressed as G1 )

G2 )

{

dV1 V1,min(t) + V1,max(t) ) K1V1 ) K1 + dt 2 V1,max(t) - V1,min(t) V1 2

}

(10)

dV2 V2,min(t) + V2,max(t) ) K2V2 ) K2 + dt 2 V2,max(t) - V2,min(t) V2 2

}

(11)

{

where K1 and K2 are constant. The temporal evolution of the exact particle composition distribution is expressed by fV′(V1, V2, t) ) fV′(V1, V2, 0) exp[-(K1 + K2)t]

(12)

Manipulation of the moments evolution equations (6), using the growth rates (10) and (11), provides the following equation: dµi1,i2

) -µi1,i2(K1 + K2) (13) dt As in the monovariate case,2 the bivariate moments vs time relationship is an exponential law. The same exponential law is obtained computing the dimensionless moments from the exact particle composition distribution in eq 12, thus demonstrating that, in the case of linear growth rates, the results obtained by the FCMOM are exact. 3.3. Diffusion-Controlled Growth. Diffusion-controlled processes are typical in clouds, where particle radii are larger than 1 µm,3 and in crystallization, particularly in nonagitated systems. A bivariate formulation of a diffusion-controlled process is required in crystallization when dealing with particles characterized by two characteristic lengths (i.e., rodlike crystals), because the growth rates along each characteristic length depend on the characteristic lengths themselves.16,17 Diffusion-controlled growth is an application in which the moments evolution equations are not closed.2,3 In fact, the growth integrals (in the terms IG1 and IG2 in the left-hand side of eq 6) cannot be expressed directly in terms of the lower order moments. In the FCMOM, the closure is obtained using the series expansion (2).2 In this subsection, the FCMOM is applied to the growth of rodlike particles. The two particle characteristic lengths are L1 (length) and L2 (width), and an equivalent diameter Leq is defined as Leq ) [(4L1L2 + 2L2L2)/π]1/2.17 Here, an application is considered in which L1, for each particle, is larger than L2 and in which the particle growth rate is diffusion-controlled in the direction of L1 while the particle growth rate is constant in the direction of L2. When L1 is much larger than L2, Leq can be estimated by [(4L1L2)/ π]1/2. Diffusion-controlled growth rates for rodlike particles in nonagitated systems are inversely proportional to the equivalent diameter Leq.17 Therefore, the growth rate along L1 can be approximately assumed as G1 )

( )

dL1 4L1L2 ) K1 dt π

-1/2

(14)


while the growth rate along L2 is given by dL2 ) K2 (15) dt In the simulations performed, the initial dimensionless bivariate particle distribution is expressed as G2 )

fL′(L1, L2, 0) ) (1 + L1) 1(1 - L1) 1(1 + L2) 2(1 - L2) q

p

q

p2

(16)

The values of the parameters of the initial distribution function are q1 ) 2, p1 ) 4, q2 ) 3, and p2 ) 1. The internal variables’

domain is initially rectangular; the initial minimum and maximum values of the particle characteristic lengths are L1,min(0) ) 50 µm, L1,max(0) ) 100 µm and L2,min(0) ) 1 µm, L2,max(0) ) 5 µm. The constants in the growth rate laws (14) and (15) are K1 ) 1 µm2/s and K2 ) 10-3 µm/s. The simulations were stopped at the time tfin ) 500 s. The initial moments were computed from eq 16 using eq 3. The moments temporal evolution was obtained by the numerical integration of eq 6. The moving boundaries are reconstructed by solving eqs 14 and 15 applied to the boundary points.

Figure 3. Bivariate distribution of particle characteristic lengths vs particle characteristic lengths: initial condition (circles) and numerical solution at the final time by the FCMOM (solid line) for particle dissolution (constant dissolution rates, equal to K). The product between K and the final time tfin is Ktfin ) 4 × 10-1 µm. M ) 10.

Figure 4. Bivariate distribution of particle characteristic lengths vs particle characteristic lengths: numerical solution at the final time by the FCMOM (solid line) and exact solution (circles) for particle dissolution (constant dissolution rates, equal to K). The product between K and the final time tfin is Ktfin ) 4 × 10-1 µm. M ) 10.


Figure 5. Bivariate particle composition distribution vs volumes of components 1 and 2 in a particle: gamma initial condition (case 1) for particle aggregation (constant aggregation kernel).

Figure 6. Bivariate particle composition distribution vs volumes of components 1 and 2 in a particle: comparison between the numerical solution by the FCMOM of the FSE (solid line) and the exact solution of the CSE (circles) at the final time for particle aggregation (constant aggregation kernel and gamma initial condition). The final time is equal to 180 s. M ) 19.

The simulation results are plotted in Figures 1 and 2. Figure 1 shows the initial condition and the numerical solution of the bivariate distribution of the particle characteristic lengths fL′, obtained by the FCMOM using moments up to the 11th order (M ) 12). In Figure 2, the numerical solution is compared with the analytical solution at the final time, showing excellent agreement. It is noteworthy that in the final solution the shape of the internal variables’ domain is no longer rectangular.

The computational time required for the simulations is low. In fact, using a 1.6 GHz M processor (Intel Pentium) and implementing the code in Matlab, the run times are on the order of 1 s. 3.4. Particle Dissolution. In this subsection, the process of dissolution of anisotropic particles is considered. Particles dissolve and become smaller until one of the two characteristic lengths reaches a size of zero; then they disappear. Particle dissolution, even in the monovariate case, is a critical problem for moments-based methods, such as the QMOM.


Figure 7. Bivariate particle composition distribution vs volumes of components 1 and 2 in a particle: mixed gamma and exponential initial condition (case 2) for particle aggregation (constant aggregation kernel).

Figure 8. Bivariate particle composition distribution vs volumes of components 1 and 2 in a particle: comparison between the numerical solution by the FCMOM of the FSE (solid line) and the exact solution of the CSE (circles) at the final time for particle aggregation (constant aggregation kernel and mixed gamma and exponential initial condition). The final time is equal to 180 s. M ) 19.

However, Strumendo and Arastoopour2 demonstrated that the FCMOM provides a numerical solution in excellent agreement with the exact solution. This capability of the FCMOM is confirmed when the FCMOM is applied to the bivariate case. In the simulations performed, the initial condition of the distribution of the particle characteristic lengths is again given by eq 16, with q1 ) 2, p1 ) 4, q2 ) 3, and p2 ) 1. In addition, L1,min(0) ) 0 µm, L1,max(0) ) 2 µm, L2,min(0) ) 0 µm, and L2,max(0) ) 0.5 µm. The two dissolution rates along the two characteristic lengths are assumed to be same and equal to a constant K. The product between K and the final time tfin is Ktfin ) 4 × 10-1 µm.

The results of the numerical simulations are shown in Figures 3 and 4. In Figure 3, the initial distribution of the particle characteristic lengths and the final distribution obtained by the FCMOM are plotted; in Figure 4, the comparison between the numerical solution and the exact solution is shown. As in the monovariate case,2 the numerical solution obtained by the FCMOM is in excellent agreement with the analytical solution. The numerical solution is shown for M ) 10. The computational time required to solve the case of particle dissolution is of the same order of magnitude as in the case of diffusion-controlled growth.


Figure 9. Bivariate particle composition distribution vs volumes of components 1 and 2 in a particle: gamma initial condition for simultaneous particle aggregation (constant aggregation kernel) and growth (linear).

Figure 10. Bivariate particle composition distribution vs volumes of components 1 and 2 in a particle: comparison between the numerical solution by the FCMOM of the FSE (solid line) and the exact solution (circles) at the final time for particle aggregation (constant aggregation kernel) and growth (linear) with a gamma initial condition. The final time is equal to 50 s. M ) 15.

4. Applications to Particle Aggregation In this section, the FCMOM is applied to particle aggregation. It is assumed that each particle consists of two chemical species (bicomponent system) and the volume occupied by the two components in a particle is V1 and V2, respectively. The internal variables ξ1, ξ2 are identified with V1 and V2, respectively; therefore, V1 ) ξ1 and V1 ) ξ1. Accordingly, the particle distribution f(ξ1,ξ2,t) is identified with the particle composition distribution fV(V1,V2,t).

In the monovariate case, Strumendo and Arastoopour2 applied the FCMOM to particle aggregation using different aggregation models. Specifically, both a Finite version of the Smoluchowski Equation (FSE) and the Oort-Hulst equation (OHE) were numerically solved by the FCMOM. A requirement for the application of the FCMOM is that the aggregation model considered must be defined in a finite domain (as the FSE and the OHE models). Similarly, in the bivariate case, aggregation models satisfying the requirement that they are defined in a finite


domain can be solved by the bivariate formulation of the FCMOM. Specifically, in this subsection, a bivariate FSE is defined and solved by the FCMOM. The simulation results are compared with analytical solutions of the bivariate classical Smoluchowski equation (CSE)19,20 for the case of constant kernel. The CSE cannot be solved directly using the FCMOM, because the CSE is defined in an infinite domain and in the CSE perturbations spread with infinite speed.21 The bivariate formulation of the CSE is ∂fV(V1, V2, t) 1 V1 V2 ) 0 0 β(V1 - η1, η1, V2 - η2, η2) × ∂t 2 fV(V1 - η1, V2 - η2, t) fV(η1, η2, t) dη1 dη2 -

∫ ∫

fV(V1, V2, t)

∫∫ ∞

0

∞

0

β(V1, η1, V2, η2) fV(η1, η2, t) dη1 dη2 (17)

In eq 17, the first term on the right-hand side represents the birth of particles due to aggregation; the second term on the right-hand side represents the death of particles due to aggregation. β is the aggregation kernel. The bivariate FSE generalizes the monovariate FSE.2 The main difference of the FSE model, when compared to the CSE model, is that in the FSE model aggregations leading to particles larger than V1,max or V2,max are neglected. Therefore, in the FSE model, the boundaries V1,min, V1,max, V2,min, and V2,max do not move with respect to time and remain equal to the values prescribed initially (at time zero). The FSE model may be expressed as follows: ∂fV(V1, V2, t) 1 ) H(V1 - 2V1,min) H(V2 - 2V2,min) × ∂t 2

∫

V1-V1,min

V1,min

∫

V2-V2,min

V2,min

β(V1 - η1, η1, V2 - η2, η2) ×

fV(V1 - η1, V2 - η2, t) fV(η1, η2, t) dη1 dη2 fV(V1, V2, t) H[(V1,max - V1,min) - V1] H[(V2,max - V2,min) - V2] ×

∫

V1,max-V1

V1,min

∫

V2,max-V2

V2,min

β(V1, η1, V2, η2) fV(η1, η2, t) dη1 dη2 (18)

In eq 18, the first term on the right-hand side represents the birth of particles due to aggregation; the second term on the right-hand side represents the death of particles due to aggregation. H is the Heaviside step function. The FSE model (18) is a specific case of the general bivariate PBE (5), in which the terms of birth and death (right-hand side of eq 5) are given by the right-hand side of eq 18. Therefore, when the FCMOM is applied to the FSE model (18), the moments evolution equations to be solved are eq 6, in which, on the left-hand side, the only term considered is the accumulation term, and in the right-hand side the FSE model (18) is introduced in the moments transformation. Closure is again obtained using the series expansion (2). Simulations were performed with a constant kernel β ) β0, and two cases were considered, using two different initial conditions. In the first case, the initial condition of the particle composition distribution fV(V1,V2,t) was a gamma distribution (case 2 in Gelbard and Seinfeld19):

(

fV(V1, V2, 0) ) 16

)( (

)

Nin V1 V2 × V1,av,inV2,av,in V1,av,in V2,av,in V1 V2 exp -2 (19) exp -2 V1,av,in V2,av,in

) (

)

In the second case, the initial condition of the particle composition distribution fV(V1,V2,t) was a mixed gamma and exponential distribution (case 3 in Gelbard and Seinfeld19):

(

fV(V1, V2, 0) ) 4

)( ) (

)

V1 Nin V1 exp -2 × V1,av,inV2,av,in V1,av,in V1,av,in V2 (20) exp V2,av,in

(

)

In eqs 19 and 20, V1,av,in and V2,av,in are the initial average values of V1 and V2 (average with respect to the initial distribution function in the [0,∞] domain), both equal to 4.189 × 10-15 m3, and Nin, the initial number of particles per unit volume, is 109/ 4.189 particles/m3.22 In the simulations, the kernel β0 ) 1.8 × 10-10 m3/s.22 The initial moments were computed from eqs 19 and 20 using eq 3; the scale factor is fsc )

Nin V1,av,inV2,av,in

Figures 5 and 6 and Figures 7 and 8 show the simulation results for cases 1 and 2, respectively. In both cases, the minimum values of each component volume were set equal to zero. As far as case 1, Figure 5 shows the initial condition of the particle composition distribution; the ratio V1,max/V1,av,in ) V2,max /V2,av,in ) 10. Figure 6 shows the comparison, at the final time, between the numerical solution by the FCMOM and the analytical solution;19 the numerical solution was obtained using moments up to the 18th order (M ) 19). As far as case 2, Figure 7 shows the initial condition of the particle composition distribution, in which the ratio V1,max/V1,av,in ) V2,max/V2,av,in ) 10. Figure 8 shows the comparison, at the final time, between the numerical solution by the FCMOM (M ) 19) and the analytical solution.19 In both cases, the numerical solution obtained by the FCMOM is in excellent agreement with the analytical solution at the final time, when the ratio N/Nin is approximately 0.2. As expected, when the ratios V1,max/V1,av,in and V2,max/V2,av,in were increased, the numerical solution of the FSE converged to the analytical solution of the CSE. In all the simulations performed for particle aggregation, the computational time required was a few seconds. 5. Application to Simultaneous Particle Aggregation and Growth In this section, the FCMOM is applied to a bivariate process, in which both aggregation and growth are present simultaneously. More specifically, the case with constant aggregation kernel and with linear growth rate laws is considered, and the model used to define the aggregation process is the bivariate formulation of the FSE (18). As in section 4, also in this section each particle consists of two components and the internal variables are the volumes, V1 and V2, of each component in a particle. Again, the particle distribution f(ξ1,ξ2,t) is identified with the particle composition distribution fV(V1,V2,t). In the application of the FCMOM to the simultaneous aggregation and growth process, the equations to be solved are eq 6 and the moving boundary conditions. In eq 6, an aggregation model (here, the FSE model) is used to define the birth and death of particles due to aggregation (right-hand side of eq 6). Differently from the solution of the FSE model for pure aggregation (section 4), in simultaneous aggregation and growth the minimum and maximum volumes are to be considered functions of time. The moving boundary conditions are obtained applying the growth rate laws (10) and (11) to the minimum and maximum volumes of each component.


In the simulations performed, the initial particle composition distribution was eq 19. The values of V1,av,in, V2,av,in, Nin, β0, and fsc are the same as used in section 4. In the linear growth rates G1 ) dV1/dt ) K1V1 and G2 ) dV2/dt ) K2V2, the constants assume the following values: K1 ) 6/21600 s-1 and K2 ) 600/ 21600 s-1. Correspondingly, the dimensionless parameters λ1 ) K1/(β0Nin) and λ2 ) K2/(β0Nin), which are the ratios between the characteristic time of growth of each component vs the characteristic time of aggregation, are λ1 ) 0.006 and λ2 ) 0.646, respectively. Figures 9 and 10 show the simulation results. Particularly, in Figure 9 the initial condition is plotted; the initial maximum volumes were set such that V1,max,in/V1,av,in ) V2,max,in/V2,av,in ) 6. Figure 10 shows the numerical solution obtained by the FCMOM and the corresponding analytical solution.19 Once again, the solution provided by the FCMOM (obtained with M ) 15) represents very nicely the exact solution. In the case of simultaneous particle aggregation (constant kernel) and growth (linear), the computational time required to perform the simulations was approximately 11 s. 6. Conclusions The FCMOM, developed by Strumendo and Arastoopour2 to solve PBE and previously validated for monovariate PBE, was extended in the present study for solving bivariate PBE. The FCMOM was applied to particle growth (constant, linear, diffusion controlled), to particle dissolution, to particle aggregation, and to simultaneous particle growth and aggregation. In all the cases considered, the solution converged to the exact solution. An important advantage of the FCMOM is that, even in the bivariate applications, it is capable of reconstructing the bivariate particle distribution function from the moments. The reconstruction is accurate, obtained through a simple relationship and valid in both transient and steady-state conditions. Additionally, the FCMOM maintains in the bivariate applications the computational efficiency demonstrated in the monovariate case;2 this feature is obviously of primary importance from the perspective of more complex applications in which the PBE are coupled with CFD (computational fluid dynamics). Finally, in the FCMOM, the internal variables always assume finite values and the domains of the internal variables are clearly defined within boundaries, whose evolution can be tracked by moving boundary conditions; this assures that the values assumed by the internal variables are always realistic. Nomenclature B ) term of birth of particles c ) coefficients in the bivariate particle distribution function series expansion D ) term of death of particles f ) bivariate particle distribution function fL ) distribution function of the particle characteristic lengths fV ) particle composition distribution function G ) particle growth/dissolution rate H ) Heaviside step function K ) constant in the growth/dissolution rate law L ) particle characteristic length Leq ) particle equivalent diameter M ) number of terms in the truncated series expansion approximating the particle bivariate distribution (equal to the order of the moments used in a simulation) m ) summation index

n ) summation index N ) number of particles per unit volume p ) parameter in the initial bivariate distribution function q ) parameter in the initial bivariate distribution function t ) time V ) volume Greek Symbols β ) aggregation kernel β0 ) constant aggregation kernel η ) volume λ ) dimensionless parameter (ratio between the characteristic time of growth and the characteristic time of aggregation) µ ) moment of the bivariate distribution function ν ) summation index σ ) summation index ξ ) internal variable φ ) Legendre trial functions Subscripts/Superscripts/Accents av ) average value fin ) final value i ) moment order in ) initial value min ) minimum value max ) maximum value sc ) scale factor - ) dimensionless quantities ′ ) denotes function of the dimensionless internal variable ξj

Literature Cited (1) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Systems in Engineering; Academic Press: London, 2000. (2) Strumendo, M.; Arastoopour, H. Solution of PBE by MOM in Finite Size Domains. Chem. Eng. Sci. 2008, 63 (10), 2624. (3) McGraw, R. Description of Aerosol Dynamics by the Quadrature Method of Moments. Aerosol Sci. Technol. 1997, 27, 255. (4) Barrett, J. C.; Webb, N. A. A Comparison of Some Approximate Methods for Solving the Aerosol General Dynamic Equation. J. Aerosol Sci. 1998, 29, 31. (5) Marchisio, D. L.; Pikturna, J. T.; Fox, R. O.; Vigil, R. D. Quadrature Method of Moments for Population-Balance Equations. AIChE J. 2003, 49, 1266. (6) Marchisio, D. L.; Vigil, R. D.; Fox, R. O. Quadrature Method of Moments for Aggregation-Breakage Processes. J. Colloid Interface Sci. 2003, 258, 322. (7) Wright, D. L.; McGraw, R.; Rosner, D. E. Bivariate Extension of the Quadrature Method of Moments for Modeling Simultaneous Coagulation and Sintering of Particle Populations. J. Colloid Interface Sci. 2001, 236, 242. (8) McGraw, R.; Wright, D. L. Chemically Resolved Aerosol Dynamics for Internal Mixtures by the Quadrature Method of Moments. J. Aerosol Sci. 2003, 34, 189. (9) Yoon, C.; McGraw, R. Representation of Generally Mixed Multivariate Aerosols by the Quadrature Method of Moments: I. Statistical foundation. J. Aerosol Sci. 2004, 35, 561. (10) Yoon, C.; McGraw, R. Representation of Generally Mixed Multivariate Aerosols by the Quadrature Method of Moments: II. Aerosol dynamics. J. Aerosol Sci. 2004, 35, 577. (11) Diemer, R. B.; Olson, J. H. A Moment Methodology for Coagulation and Breakage Problems: Part 1-Analytical Solution of the Steady State Population Balance. Chem. Eng. Sci. 2002, 57, 2193. (12) Diemer, R. B.; Olson, J. H. A Moment Methodology for Coagulation and Breakage Problems: Part 2-Moments Models and Distribution Reconstruction. Chem. Eng. Sci. 2002, 57, 2211. (13) Diemer, R. B., Jr.; Olson, J. H. Bivariate Moment Methods for Simultaneous Coagulation, Coalescence and Breakup. J. Aerosol Sci. 2006, 37, 363. (14) Marchisio, D. L.; Fox, R. O. Solution of Population Balance Equations Using the Direct Quadrature Method of Moments. J. Aerosol Sci. 2005, 36, 43.

Ind. Eng. Chem. Res., Vol. 48, No. 1, 2009 273 (15) Hulburt, H. M.; Katz, S. Some Problems in Particle Technology. Chem. Eng. Sci. 1964, 19, 555. (16) Puel, F.; Fevotte, G.; Klein, J. P. Simulation and Analysis of Industrial Crystallization Processes through Multidimensional Population Balance Equations. Part 2: a Study of Semi-Batch Crystallization. Chem. Eng. Sci. 2003, 58, 3729. (17) Puel, F.; Fevotte, G.; Klein, J. P. Simulation and Analysis of Industrial Crystallization Processes through Multidimensional Population Balance Equations. Part 1: a Resolution Algorithm Based on the Method of Classes. Chem. Eng. Sci. 2003, 58, 3715. (18) McMurry, P. H.; Wilson, J. C. Droplet Phase (Heterogeneous) and Gas Phase (Homogeneous) Contributions to Secondary Ambient Aerosol Formation as Functions of Relative Humidity. J. Geophys. Res. 1983, 88, 5101.

(19) Gelbard, F. M.; Seinfeld, J. H. Coagulation and Growth of a Multicomponent Aerosol. J. Colloid Interface Sci. 1978, 63, 472. (20) Smoluchowski, M. Z. Versuch Einer Mathematischen Theorie der Koagulationskinetik Kolloider Losungen. Z. Phys. Chem. 1917, 92, 129. (21) Dubovski, P. B. A “Triangle” of Interconnected Coagulation Models. J. Phys. A: Math. Gen. 1999, 32, 781. (22) Scott, W. T. Analytic Studies of Cloud Droplet Coalescence I. J. Atmos. Sci. 1968, 25, 54.

ReceiVed for reView February 17, 2008 ReVised manuscript receiVed August 25, 2008 Accepted August 28, 2008 IE800272A

Solution of Bivariate Population Balance Equations Using the Finite

Recommend Documents