Bivariate direct quadrature method of moments for coagulation and sintering of particle populations

Aerosol Science 37 (2006) 1562 – 1580 www.elsevier.com/locate/jaerosci

Bivariate direct quadrature method of moments for coagulation and sintering of particle populations R.O. Fox∗ Department of Chemical and Biological Engineering, 2114 Sweeney Hall, Iowa State University, Ames, IA 50011-2230, USA Received 3 December 2005; received in revised form 6 March 2006; accepted 10 March 2006

Abstract The direct quadrature method of moments (DQMOM) is applied to a bivariate population balance equation governing the joint number density function n(v, a) of particle volume and area. Using a test case proposed by Wright et al. [2001. Bivariate extension of the quadrature method of moments for modeling simultaneous coagulation and sintering of particle populations. Journal of Colloid and Interface Science, 236, 242–251], it is shown that DQMOM reproduces the results found by using bivariate extensions of QMOM. The effects of the choice of bivariate moments and the order (N) of the quadrature approximation are also discussed. In general, it is found that for fixed N all non-singular sets of bivariate moments yield nearly identical results as reported by Wright et al. However, particular sets leading to lower condition numbers for the linear system used to compute the source terms in DQMOM are preferable for numerical reasons. For any value of N
3, the bivariate moment predictions are very similar. Overall, it is found that the coagulation and sintering kernels used in this test case are easily handled by lower-order moment methods, and are relatively insensitive to the choice of bivariate moments used for DQMOM. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Direct quadrature method of moments; Bivariate population balance; Coagulation; Sintering

1. Introduction The quadrature method of moments (QMOM) was introduced by McGraw (1997) to solve for the moments of a particle size distribution (PSD) for arbitrarily complex growth and coagulation laws. Starting with transport equations for the univariate moments, the source terms due to physical processes such as growth and coagulation are closed by expressing them in terms of the weights wn and abscissas vn . The latter are related to the k-order moment Mk by a quadrature formula: Mk =

N


wn vnk .

n=1

∗ Tel.: +1 515 294 9104.

E-mail address: [email protected]. 0021-8502/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.jaerosci.2006.03.005

(1)

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

1563

As shown by McGraw (1997) the N weights and N abscissas are uniquely determined by a set of 2N moments, and a robust algorithm exists for determining (wn , vn ) from Mk . This fact has been exploited by several groups (see, for example, Diemer & Ehrman, 2005; Marchisio, Pikturna, Fox, Vigil, & Barresi, 2003; Marchisio, Vigil, & Fox, 2003) to solve univariate PSDs undergoing any combination of nucleation, growth, aggregation and breakage. More recently, Marchisio and Fox (2005) have shown that, instead of solving transport equation for the moments, it is also possible to derive and solve transport equations for the weights and abscissas directly. This alternative method is known as direct QMOM or DQMOM. For univariate PSD it is straightforward to show that QMOM and DQMOM (with integer moments) are mathematically equivalent (Marchisio & Fox, 2005). However, for the univariate case DQMOM can also be used with independent, noninteger moments, which may offer advantages for certain applications (Diemer & Ehrman, 2005; Upadhyay & Ezekoye, 2006). Independently, Piskunov and Golubev (1999) developed a closely related method (called the GA method in Piskunov & Golubev, 2002) that solves ordinary differential equations directly for variables that are equivalent to the weights and abscissas used in quadrature methods. In principle, the extension of QMOM to a multivariate PSD is straightforward. For example, the moments of a bivariate PSD are related to the weights and abscissas (e.g., volume and area) by

Mkl =

N


wn vnk anl .

(2)

n=1

As discussed by Wright, McGraw, and Rosner (2001), the main difficulty in applying QMOM to bivariate problems is the lack of an efficient numerical algorithm for determining the weights and abscissas from the moments. Despite this difficultly, these authors have shown that bivariate QMOM yields accurate results for simultaneous coagulation and sintering in spatially homogeneous systems at a fraction of the cost of methods that solve for the entire distribution function. However, the use of QMOM for spatially inhomogeneous systems still relies on solving transport equations for the bivariate moments. Thus, at every time step and in every grid cell in the computational domain, Eq. (2) must be inverted to find the weights and abscissas. This inversion procedure is relatively expensive from a computational standpoint. Moreover even for the univariate case, due to discretization errors when solving moment-based transport equations numerically, the moments can fail to be realizable (e.g., when M2 M12 ), making quadrature inversion impossible. In contrast, DQMOM solves spatial transport equations for the weights and abscissas directly and, as long as the weights remain non-negative in the numerical algorithm, the moments will always be realizable. Thus, DQMOM offers a great potential for lowering the computational cost in computational fluid dynamic (CFD) simulations of evolving particle populations (see Fan, Marchisio, & Fox, 2004; Zucca, Marchisio, Barresi, & Fox, 2006 for examples). Piskunov and Golubev (2002) have also proposed an extension of the GA method to bivariate coagulation kinetics, which presumably should yield ordinary differential equations equivalent to DQMOM for the spatially homogeneous case (this is not shown, however, in Piskunov & Golubev, 2002). Thus, the two principal advantages of using the DQMOM methodology are that (1) it is easily applied to spatially inhomogeneous transport equations wherein new terms arise due to spatial gradients (Marchisio & Fox, 2005), size-dependent convection (Fan et al., 2004), and size-dependent diffusion (Upadhyay & Ezekoye, 2006); and (2) it keeps all of the moments realizable when used in discretized transport models (e.g., for atmospheric and climate modeling). The goal of the present work is to demonstrate that DQMOM is equivalent to QMOM for modeling simultaneous coagulation and sintering of spatially homogeneous bivariate particle populations. For this purpose, we will use the example problem of coagulation and sintering proposed by Wright et al. (2001). We will also investigate how the choice of the bivariate moments affects the weights and abscissas and the numerical properties of DQMOM. As discussed in the recent work of Diemer and Olson (2006), efficient solvers based on moment methods will open the door to modeling spatially inhomogeneous bivariate processes, which are still too costly to treat using sectional methods coupled with CFD. The remainder of this work is arranged as follows. In Section 2 we introduce the transport equation for the bivariate PSD, and the coagulation and sintering rates used by Wright et al. (2001). In Section 3 we review the DQMOM transport equations for the bivariate case and discuss the linear system used to solve for the source terms for the weights and abscissas. In Section 4 we present and discuss the results for coagulation and sintering. Conclusions are drawn in Section 5. Details on the DQMOM derivation are included in the Appendix.

1564

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

2. Bivariate population balance equation A transport equation for the joint volume, area number density function n(v, a; x, t) can be expressed as jt n + u · jx n − jx · (Djx n) + jv (Rv n) + ja (Ra n) = Q,

(3)

where u is a solenoidal convection velocity (jx · u = 0), D is a diffusion coefficient, Rv is the volume-growth rate, Ra is the area-growth rate, and Q is the coagulation term. The form of the temporal and spatial transport terms will influence the source terms in DQMOM. Here we include convection and an isotropic diffusion term in order to illustrate their effects. Nevertheless, DQMOM can be applied to more general transport equations. The extension of DQMOM to treat size-dependent transport processes is straightforward (Fan et al., 2004; Upadhyay & Ezekoye, 2006), but will not be needed in this work. Using standard assumptions (Marchisio & Fox, 2005; Wright et al., 2001), we can write Q = Q+ + Q− where   1 ∞ ∞ Q+ =
(v − v ∗ , v ∗ , a − a ∗ , a ∗ )n(v − v ∗ , a − a ∗ )n(v ∗ , a ∗ ) dv ∗ da ∗ , (4) 2 0 0  ∞ ∞ −
(v, v ∗ , a, a ∗ )n(v, a)n(v ∗ , a ∗ ) dv ∗ da ∗ . (5) Q =− 0

0

The coagulation rate
, and the growth rates Ra and Rv can generally be functions of both volume and area, and DQMOM can easily treat such cases. However, following Wright et al. (2001), we will assume here that
depends only on volume:
(v, u) = K(v −1/Df + u−1/Df )(v 1/Df + u1/Df ),

(6)

where K is a proportionality constant and Df is the particle fractal dimension. (For Brownian coagulation of spherical particles, Df = 3.) Likewise, following Wright et al. (2001), we will take Rv = 0 and use a linear expression (in a) for Ra (Koch & Friedlander, 1990): Ra (v, a) =

1 [amin (v) − a] tf

(7)

where amin (v) = (6v/)2/3 and tf is the characteristic sintering time. Nevertheless, one of the key advantages of DQMOM (and quadrature methods in general) is that arbitrarily complex nonlinear growth rates can be treated without any additional complications. 3. DQMOM approximation of population balance equation DQMOM approximates the density function by weighted delta functions in v-a phase space (Fox, 2003; Marchisio & Fox, 2005): n(v, a) =

N


wn (v − vn )(a − an ).

(8)

n=1

Note that in this formulation, the weights wn and abscissas (vn , an ) are fields (i.e., they depend on space as well as time). The DQMOM approximation for the bivariate moments at any time instant and spatial location is found directly from Eq. (8): ∞ ∞

 Mkl ≡ 0

0

v k a l n(v, a) dv da =

N


wn vnk anl .

(9)

n=1

The fundamental idea behind QMOM is that we should determine the weights and abscissas such that as many bivariate moments computed from Eq. (9) as possible agree with same moments computed from the transport equations (found from Eq. (3)): jt Mkl + u · jx Mkl − jx · (Djx Mkl ) = S(k, l),

(10)

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

1565

where S(k, l) is the moment source term (see Eq. (48) in Appendix for the exact form). Note that there are a total of N weights and 2N abscissas in Eq. (9) and (equivalently) 3N moment source terms in Eq. (10). We will thus need to choose 3N independent bivariate moments to determine the source terms in the DQMOM transport equations for the weights and abscissas. We will return to the subject of how to choose the moments in Section 3.2. 3.1. DQMOM transport equations Application of DQMOM (see Appendix) results in closed transport equations for the number density, volume density, and area density, respectively, of the form: jt wn + u · jx wn = jx · (Djx wn ) + An ,

(11)

jt (wn vn ) + u · jx (wn vn ) = jx · (Djx wn vn ) + Bn ,

(12)

and jt (wn an ) + u · jx (wn an ) = jx · (Djx wn an ) +

wn [amin (vn ) − an ] + Cn , tf

(13)

where An , Bn , and Cn are source terms that are found by forcing Eqs. (11)–(13) to agree with Eqs. (9) and (10). These source terms are computed by solving a linear system: N


l k l−1 [(1 − k − l)vnk anl An + kv k−1 n an Bn + lv n an Cn ] = P (k, l) + D(k, l),

(14)

n=1

where the term on the right-hand side due to coagulation is given by P (k, l) =

N N 1

wn wq [(vn + vq )k (an + aq )l − vnk anl − vqk aql ]
(vn , vq ), 2

(15)

n=1 q=1

and the term due to diffusion by D(k, l) =

N


l−1 k l−2 [k(k − 1)vnk−2 anl wn Dvvn + 2klv k−1 n an wn Dvan + l(l − 1)vn an wn Daan ],

(16)

n=1

wherein the gradient-dependent terms are defined by Dvvn = D(jx vn ) · (jx vn ),

Dvan = D(jx vn ) · (jx an ),

and

Daan = D(jx an ) · (jx an ).

(17)

As described below, the linear system in Eq. (14) will be well defined for any independent set of 3N bivariate moments. When solving the linear system in Eq. (14) numerically, it is useful to note the scaling properties of the coefficient matrix on the left-hand side. For coagulation problems, the abscissas can grow to be very large and it is helpful to rescale the problem to avoid round-off error. If we introduce the scaling factors vs and as defined by vs = max vn n

and

as = max an , n

(18)

the scaled abscissas can be defined by vn∗ = vn /vs and an∗ = an /as . We can then divide Eq. (14) by vsk asl , and find that it has exactly the same form if we define Bn = vs Bn∗

and

Cn = as Cn∗ ,

(19)

where Bn∗ and Cn∗ are found by solving Eq. (14) using the scaled abscissas in place of the unscaled values. From a practical standpoint, it is also important to note that the left-hand side of Eq. (14) does not involve the weights wn . This implies that the condition number of the linear system will be independent of the particle number density (in fact, the system will be well defined for wn = 0). This property is especially important when dealing with particulate systems that are generated by nucleation or condensation (Zucca et al., 2006).

1566

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

Before continuing, it is interesting to remark on the properties of Eqs. (11)–(13). First, we see that the convection terms from Eq. (3) (i.e., the terms involving u, Rv and Ra ) appear in a closed, uncoupled form in the transport equations for the weights and abscissas. Thus, in the absence of diffusion and coagulation, the weights remain unchanged along streamlines following the local fluid velocity and the abscissas evolve according to the “local” values of Rv and Ra . Readers familiar with Lagrangian particle methods will recognize this property as one also shared by particle methods (Fox, 2003). On the other hand, when diffusion and coagulation are included, the abscissas and weights are coupled through the “non-local” source terms. For the coagulation terms, non-local interactions follow our intuition of particleparticle interactions (albeit here represented by abscissa–abscissa coupling). However, the diffusion terms might at first appear to be unexpected. The origin of these terms can be understood by considering the effect of random particle motion along the direction of a mean gradient in (say) the volume. If initially the particles are aligned spatially according to their volume, by definition the variance of the volume will be null at every location but the mean volume will vary with location. Introducing random diffusive steps will then yield an increase in the local variance. It is this process that generates the source term D(k, l) appearing in Eq. (14). The derivation procedure used in DQMOM provides the correct form for the diffusive source terms directly. If these terms were neglected (e.g., if DQMOM source terms were derived only by looking at the homogeneous equations), then the bivariate moments would be predicted incorrectly. For cases where the convection and/or diffusion terms are size dependent (Fan et al., 2004; Upadhyay & Ezekoye, 2006), it is even more crucial to systematically apply the DQMOM derivation procedure to find all of the necessary source terms. Eqs. (11)–(13) can be solved with appropriate initial and boundary conditions to find the fields wn (x, t), vn (x, t), and an (x, t) appearing in Eq. (8). An example of how this is done is given in Wang and Fox (2004). Note that once initial and boundary conditions have been established for the weights and abscissas, the latter are determined at all future times and locations by Eqs. (11)–(13). This is the principal difference between DQMOM and QMOM, where in the latter Eq. (10) is solved for the bivariate moments, and a numerical algorithm is needed to convert the updated moments to weights and abscissas. Similar to DQMOM, another “direct” method (Jacobian matrix transformation) has been proposed by McGraw and Wright (2003). For the homogeneous case, JMT and DQMOM are mathematically equivalent. An important difference, however, is the term D(k, l), which appears in DQMOM, but not in JMT. As noted above, it is this term that correctly accounts for the effects of spatial gradients in the weights and abscissas for inhomogeneous cases, and avoids the need to solve the transport equations for the moments (Eq. (10)) directly. 3.2. DQMOM linear system The DQMOM approximation of the bivariate population balance equation is given by the transport equations for the weights and abscissas (Eqs. (11)–(13)). The source terms for these equations are found by solving the linear system given in Eq. (14). The exact form of the linear system depends on the choice of bivariate moments. This choice, in turn, will determine whether the system is well defined in the sense that the coefficient matrix is non-singular. Even for cases where the matrix is non-singular, when N is larger than 5 or 6 the matrix can be poorly conditioned (even for what would otherwise appear to be “good” choices of moments). In practice, this is usually due to round-off error and it can lead to artificial “stiffness” in the set of differential equations when in reality the equations are well behaved. Based on our experience using DQMOM, we have found that the round-off error can be completely eliminated by using the iterative improvement algorithm described in Section 2.5 of Press, Teukolsky, Vetterling, and Flannery (1992). Normally, the solution to Eq. (14) can be found to machine precision using at most three iterations. In all of the results reported in this work, we have used the singular value decomposition (SVD) method described in Press et al. (1992) to solve the linear system. Use of this method allows us to track easily the condition number of the matrix with no additional overhead. As with the lower–upper (LU) decomposition method, the SVD method is very efficient when combined with the iterative improvement algorithm, so very little additional expense is required to get accurate solutions to Eq. (14). Many choices of 3N moments can be shown to yield a non-singular coefficient matrix in Eq. (14). Strategies for choosing moments have been discussed by Wright et al. (2001) and Upadhyay and Ezekoye (2006). In this work, results for the different moment sets given in Tables 1–3 are reported. In general, we find that the choice of the moments (for fixed N > 2) does not strongly affect the accuracy of the moment predictions for this problem. It should be noted, however, that some choices yield better-conditioned coefficient matrices than others. The condition number of a matrix is defined as the ratio of the largest and smallest singular values (Press et al., 1992). All else being equal, it is thus preferable to choose moments that yield small condition numbers. For example, Set A4 in Table 3 has a condition

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

1567

Table 1 Set of six moments employed with N = 2 Set A2 k

l

1

0

0

2

0

3

1 3 2 3

4

1

0

5

0

6

0

1 3 2 3

0

Table 2 Set of nine moments employed with N = 3 Set A3 k

l

1

0

0

2

0

3

1 3 2 3

4

1

0

5

0

6

4 3 5 3

7

0

8

0

1 3 2 3

9

0

1

0

0

Table 3 Six sets of twelve moments employed with N = 4 Set A4

Set B4

Set C4

Set D4

Set E4

Set F4

k

l

k

l

k

l

k

l

k

l

k

l

1

0

0

0

0

0

0

0

0

0

0

0

0

2

0

1 3 2 3

1 3 2 3

0

0

0

1 3 2 3

4

1

0

1

0

1

1

1

0

0

1

1 4 1 3 1 2

0

0

1 3 2 3

0

0

1 3 2 3

0

3

1 3 2 3

5

0

4 3 5 3

4 3 5 3

0

0

0

0

4 3 5 3

1

0

4 3 5 3

0

0

4 3 5 3

0

6

4 3 5 3

2

0

7

2

0

2

0

2

2

2

0

0

2

3

0

8

7 3

0

7 3

0

0

7 3

4

0

0

0

0

1

1

1 3 2 3

0

0

0

5 3 7 3

4 3 5 3

1 3 1 3 1 3 1 3

0

0

7 3 1 3 2 3

0

0

7 3 2 3

7 3

1 3

1

0

0

1 4 1 3 1 2

4 3

0

0

1

10

0

1 3 2 3

11

0

1

0

4 3

9

12

0

1 4 3

2 3 4 3

2

0 0

1568

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

number near 1010 , while for Set C4 it is near 1016 . Hence, assuming that both sets predict the moments with equal accuracy, Set A4 would be prefered. Note also that Set C4 in Table 3 does not contain the moments (k, l) = (0, 1) and (1, 0); thus, the total volume will not be strictly conserved during coagulation and the total surface area will not be conserved in the absence of sintering. This fact, in addition to the poor condition number, might be another good reason to avoid Set C4. It is also interesting to note that Sets A, B and D–F in Table 3 yield a linear system that can be written in matrix form (showing only non-zero components) as ⎡ ⎤⎡ ⎤ ⎡ ⎤ M 1 M2 A N1 ⎢ ⎥⎢ ⎥ ⎢ ⎥ (20) ⎣ M 3 M4 ⎦ ⎣ B ⎦ = ⎣ N2 ⎦ , M5

V

C

N3

where M and V are N×N matrices and A, B and C are column vectors formed from the components An , Bn and Cn , respectively. N are column vectors formed from the right-hand side of Eq. (14). Note that A and B can be found first by solving a smaller system that corresponds to the univariate population balance equation for volume (i.e., l = 0):



M 1 M2 A N1 = . (21) B M3 M4 N2 C can then be found separately: VC = N3 − M5 A,

(22)

where V is a Vandermonde matrix (Press et al., 1992) formed from the area abscissas an . It is well known (Marchisio, Pikturna et al., 2003; Marchisio, Vigil et al., 2003) that Eq. (21) (which is equivalent to QMOM for volume) accurately predicts the 2N volume moments used to define M . An additional set of N area moments is controlled using Eq. (22). The condition numbers of M and V depend only on the volume and area abscissas, respectively. Thus, as long as the abscissas are distinct, the coefficient matrices will be non-singular. The matrix decomposition in Eq. (20) can be used to develop a strategy for choosing the moments in the case of pure coagulation. In the special case where an = cv n , the coagulation source term has the property P (k, l) = cl P (k + l, 0). It can then be easily shown that Cn = cB n , and that independent moments require distinct values of m = k + l. Using the decomposition in Eq. (21), 2N distinct values of k are needed to find An and Bn . Then, in order to ensure that Cn = cB n , it suffices to choose any subset of size N from the distinct values of k and let l take values from this subset when solving Eq. (22) (e.g., Set A4). As we will illustrate latter, this strategy leads to well-conditioned coefficient matrices. Another possible strategy can be developed by assuming that An = Aw n and using the zero-order moment to find A: A=

P (0, 0) . M00

(23)

With An known, Bn and Cn can be found separately from Eq. (14): N


kv k−1 n Bn = P (k, 0) + D(k, 0) +

n=1 N
n=1

N


(k − 1)vnk An

with k > 0,

(24)

n=1

la l−1 n Cn = P (0, l) + D(0, l) +

N


(l − 1)anl An

with l > 0.

(25)

n=1

Note that although this strategy does not fix as many moments as “standard” DQMOM, it does treat both volume and area in a symmetric, uncoupled manner. Moreover, the coefficient matrices in Eqs. (24) and (25) can be written as Vandermonde matrices of order N, which can be inverted using the routine described in Press et al. (1992). Recently, Yoon and McGraw (2004a, 2004b) proposed a method based on principal components analysis (PCA-QMOM) that also treats volume and area separately in a symmetric, uncoupled manner. These authors have also suggested fixing the weights and using only the lower-order moments to fix the abscissas. In the context of DQMOM, we shall see that

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

1569

assuming An = Awn restricts the choice of initial abscissas that can be used with a particular choice of initial weights, thereby diminishing the applicability of this assumption for system involving coagulation. Finally, we should note that with DQMOM it is not necessary to use the same set of moments for each physical phenomena. For example, because each new physical process will result in an additional term on the right-hand side of Eq. (14), we can split the moment source terms into sums over all processes: An = AP n + ADn ,

(26)

Bn = BP n + BDn ,

(27)

Cn = CP n + CDn ,

(28)

where the terms with a subscript P are found from Eq. (14) by considering only coagulation (i.e., D(k, l) = 0), etc. In our earlier work (Marchisio & Fox, 2005), we have shown that using Eqs. (23)–(25) with k, l ∈ (1, 2, . . . , N ) to find ADn , BDn and CDn for the diffusion term D(k, l) leads to realizable abscissas. On the other hand, in Section 4 we show that Eq. (20) with fractional moments works best for treating coagulation. The “combined-moment” strategy might thus be useful for treating spatially inhomogeneous cases. As shown in the Appendix, all moment sets yield the exact same source terms for sintering in the DQMOM formalism. It is thus only necessary to look at the effect of the choice of moments on the coagulation and spatial-transport terms when using DQMOM. 4. Results and discussion In order to test the numerical behavior of the DQMOM approximation, we will consider spatially homogeneous cases introduced by Wright et al. (2001) for which the weights and abscissas depend only on t and for which D(k, l) = 0. As discussed elsewhere (Marchisio & Fox, 2005), DQMOM and QMOM are equivalent for univariate and bivariate homogeneous cases when the same set of moments is used. In Wright et al. (2001) detailed comparisons are made between a selected set of 36 bivariate moments found by QMOM and those found by solving directly for n(v, a) using a bivariate sectional method. These authors consider both N = 3 and N = 12 quadrature nodes. For N = 3, different choices for the 3N = 9 moments are considered, but the choices are limited by the requirement that the moment equations must be easily invertible to find the weights and abscissas (see the original paper for details). With DQMOM, this consideration is not a factor and any linearly independent choice of bivariate moments can be used to compute the weights and abscissas. For N = 12, all 36 moments are required to determine the weights and abscissas. In Wright et al. (2001), non-linear regression is used to compute the weights and abscissas given a set of 36 moments. As discussed by these authors, relatively few iterations are required to obtain convergence. Nevertheless, the computational cost for N = 3 is approximately 100 times less than with N = 12, even though the number of unknown increases only by a factor of 4. With DQMOM, no iterations (besides those used to improve the linear solution) are required since Eq. (14) provides the source terms directly for any value of N, and thus it takes only a little more than four times longer to use N = 12 than N = 3. Nevertheless, as discussed in Section 3.2, when using N
5 one must be careful to control round-off error when solving the linear system. One of the principal conclusions drawn by Wright et al. (2001) is that QMOM with N = 3 provides good agreement with the sectional method for both the dynamic and asymptotic moments. Because DQMOM and QMOM are mathematically equivalent (if the same set of bivariate moments are used), it comes as no surprise that DQMOM also provides equally good agreement for the moments. Thus, the main focus in this work will not be on comparison with the sectional method, but rather on how the choice of bivariate moments used in DQMOM affects (1) the moment predictions and (2) the numerical properties of the linear system. 4.1. Simulation conditions The initial values for the weights and abscissas are given in Table 4, and correspond approximately to the initial moments used by Wright et al. (2001). Note that initially only one weight (w1 ) is non-zero, and the corresponding abscissas yield M10 /M00 = 1 and M01 /M00 = 25.8. The values of the other abscissas have been chosen arbitrarily and have no effect on the initial values of the moments. They do effect the values of the source terms for small times (0 t); however, the source terms quickly adjust the abscissas so that the predicted moments are independent of the

1570

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

Table 4 Initial conditions for weights and abscissas n

wn

vn

an

1 2 3 4

1 0 0 0

1 2 5 10

25.8 50 125 250

initial values of the abscissa with zero weights. This is, in fact, a significant advantage when using DQMOM since it is often difficult to find initial values for the weights and abscissas for arbitrary initial bivariate distribution functions. We should note, however, that the stiffness of the problem (as reflected in part by the condition number of the coefficient matrix) is dependent on the choice of the initial abscissas. Values that are far from “optimal” cause the system to be very stiff during a short initial period wherein they readjust. We will show examples of this behavior below. Following Wright et al. (2001), we set Df = 3 and K = 1 in Eq. (6), so that KM 00 (0) = 1. However, as noted by these authors, the results depend only weakly on the value of Df . The homogeneous system with u and D equal to zero in Eqs. (11)–(13) is solved using the initial value solver LSODA (Hindmarsh, 1983). The singular value decomposition (SVD) routines from Press et al. (1992) are used to invert the coefficient matrix in Eq. (14). Although other linear equation solvers would perform equally well, the SVD routines were chosen in order to follow the evolution of the condition number of the linear system. For all values of N that were investigated, the computation times on a standard PC are negligible, and thus will not be reported here. 4.2. Moment predictions In Fig. 1 we show results with N = 4 for the 36 bivariate moments used in Wright et al. (2001), where the bivariate moments correspond to all possible unique combinations of (k, l) ∈ (0, 13 , 23 , 1, 43 , 53 ). These plots can be compared with Fig. 2 in Wright et al. (2001) for l = 0, 13 , 53 found using QMOM with selected sets of moments. In general, the agreement between DQMOM and QMOM is, as expected, excellent. More interesting, perhaps, is the fact that all six moment sets in Table 3 yield nearly identical results. Although we do not show the results here, it was found that all linearly independent sets of bivariate moments for N
4 that we tested yield essentially identical predictions for the 36 moments. This finding suggests that, at least for this problem, when N is sufficiently large all moment choices converge to the same results. If this were true, then the only important consideration when choosing the moments is how accurately the linear system (Eq. (14)) can be solved, which will be determined primarily by the condition number of the matrix. 4.3. Condition numbers In Fig. 2 we plot the evolution of the condition number for the moment sets in Table 3. From this figure we can conclude that Set A4 might be prefered due to its significantly lower condition number. (Sets A4 and E4 differ essentially by a constant factor. This is because the condition number of Set E4 depends on a scale factor proportional to an /vn .) However, all sets except Set C4 yield acceptable condition numbers for all times. In contrast, Set C4 becomes poorly conditioned as time increases, leading to round-off errors during inversion of the matrix. Set F4 shows a similar behavior, but still yields results that are not too contaminated by numerical error. The moments in Set C4 were chosen to mimic those used by Wright et al. (2001) for N = 3. Likewise, Set F4 was formed by extending the univariate moments employed by Diemer and Ehrman (2005). Apparently these choices are not optimal from a numerical perspective. Finally, we can note from the condition-number curves that a rapid variation is observed near time zero. As discussed below, by choosing different values for the initial abscissas with zero weights, it is possible to eliminate, or at least reduce, this behavior.

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

102

101

101

100

Mkl

Mkl

1571

100 10-1

0

1

2

3

4

5 time

6

7

8

9

10-1

10

104

102

103

Mkl 101

100

2

3

4

5 6 time

7

8

9

10

5 6 time

7

8

9

10

6

7

8

9

10

102

0

1

2

3

4

5 time

6

7

8

9

101

10

105

106

104

105

Mkl

Mkl

1

Mkl

103

0

103

102

101

0

1

2

3

4

0

1

2

3

4

104

103

0

1

2

3

4

5 time

6

7

8

9

10

102

5 time

Fig. 1. DQMOM predictions with N = 4 for 36 bivariate moments for simultaneous coagulation and sintering. Each plot corresponds to a fixed value for the l-index, and the k-index increases between 0, 13 , . . . , 53 from the lowest to the highest curve. Lines correspond to moment Set A4, circles to Set B4, and squares to Set C4. Upper left: l = 0. Upper right: l = 13 . Middle left: l = 23 . Middle right: l = 1. Lower left: l = 43 . Lower right: l = 53 .

1572

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

1017

Condition number

1015

1013

1011

109 0

5

10 time

15

20

Fig. 2. Evolution of the condition number for moment sets in Table 3. Solid: Set A4. Circle: Set B4. Square: Set C4. Diamond: Set D4. Triangle up: Set E4. Triangle down: Set F4.

4.4. Weights and abscissas In Figs. 3 and 4 we have plotted the weights and abscissas for selected moment sets from Table 3. It is remarkable that all sets (only Sets A4, C4 and F4 are shown in the figures) except Sets C4 and F4 yield nearly identical weights and abscissas. In contrast, the weights and abscissas found from Set C4 (Fig. 3) (and less so with Set F4 in Fig. 4) are significantly different, even though the moment predictions (Fig. 1) are nearly identical to the other sets. From the abscissas, it can be observed that some of them (i.e., those with null weight at time zero) vary rapidly near time zero, before reaching “optimal” values lying on smoothly increasing curves. By tracing these curves back to the origin, it is possible to obtain “optimal” initial values for the abscissas that eliminate the stiffness of the system. Note that these values appear to be independent of the chosen moment set. 4.5. Dependence on N In Fig. 5 results for the 36 bivariate moments found by varying N from 2 to 4 are shown for Sets A(2–4). (Results for the dependence on N for other moment sets are consistent with Set A, and thus not shown.) Except for the highest-order moments, the results with N = 3 or 4 are nearly indistinguishable. In contrast, the results for N = 2 are only accurate for the low-order moments. Thus, as found in earlier work for univariate cases (Marchisio, Vigil et al., 2003), we can conclude that N
3 yields accurate results for the bivariate case studied by Wright et al. (2001). 4.6. Effect of other parameters Using the moment sets in Table 3, we have verified that DQMOM reproduces the QMOM results in Wright et al. (2001) for other cases (i.e., different values of Df and tf ). Since these authors have validated the (asymptotic) QMOM results with the Monte-Carlo results of Tandon and Rosner (1999), we can conclude that DQMOM yields

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

1573

100

wn

10-1

10-2

10-3

10-4 0

5

0

5

10 time

15

20

103

vn

102

101

100

10 time

15

20

Fig. 3. Weights (top) and volume abscissas (bottom) for moment sets from Table 3. Solid lines: Set A4. Dashed lines: Set C4. Circle: node 1. Square: node 2. Diamond: node 3. Triangle: node 4.

excellent results for homogeneous cases using the coagulation and sintering kernels considered in this work. Our findings presented above for the independence on moment choice for sufficiently large N, for the condition number, and for the dependence on N were also found to hold for other parameter choices.

1574

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

100

wn

10-1

10-2

10-3

10-4

0

5

0

5

10 time

15

20

104

an

103

102

101

10 time

15

20

Fig. 4. Weights (top) and area abscissas (bottom) for moment sets from Table 3. Solid lines: Set A4. Dashed lines: Set F4. Circle: node 1. Square: node 2. Diamond: node 3. Triangle: node 4.

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

102

101

101

100

Mkl

Mkl

1575

100 10-1

0

1

2

3

4

5 6 time

7

8

9

10-1

10

1

2

3

4

5

6

7

8

9

10

6

7

8

9

10

6

7

time 104

102

103

Mkl

Mkl

103

101

100

102

0

1

2

3

4

5 6 time

7

8

9

101

10

105

106

104

105

103

102

101

0

1

2

3

4

5 time

Mkl

Mkl

0

104

103

0

1

2

3

4

6 5 time

7

8

9

10

102 0

1

2

3

4

5

8

9

10

time

Fig. 5. DQMOM predictions for 36 bivariate moments for simultaneous coagulation and sintering with moment Set A(2-4). Lines correspond to N = 4, circles to N = 3, and squares to N = 2. See Fig. 1 for definition of plots and curves.

1576

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

vn

10

1

0

0.2

0.4

0.6

0.8

1

time Fig. 6. Volume abscissas for N = 3 found by assuming all weights are equal. Circle: node 1. Square: node 2. Diamond: node 3. Abscissas v1 and v2 nearly coincide at t = 0.22.

4.7. Constrained source terms We next attempted to generate results using Eqs. (23)–(25) and wn =1/N (i.e., equal weights). It was found, however, that for N
3 the condition number of the system depends strongly on the initial abscissas, and that the system becomes singular for many choices of the initial conditions. (For N = 2, the system is well behaved.) For example, we show in Fig. 6 a set of initial conditions that lead to nearly singular behavior with N =3. It can be seen in the figure that v1 and v2 approach each other near t = 0.22, before moving apart. If a smaller initial value were used for v2 , these two abscissas would coincide at a particular time, generating a singular coefficient matrix. As N increases, the problem of coinciding abscissas becomes even more difficult to avoid. The root of this problem lies in the limited range of moments that can be reproduced by abscissas with fixed weights. For N = 2, only the first two moments are controlled, and this can always be done with two abscissas regardless of the weights. However, for N
3, it is possible to have a realizable third-order moment that cannot be reproduced by three real-valued abscissas with fixed weights. As seen here for N = 3, initial conditions for which v2 is too small will lead to such a problem, and DQMOM with fixed weights will fail. This failure underscores the remarkable fact that with variable weights (i.e., “standard” DQMOM) a consistent set of real-valued abscissas can be found for any realizable set of independent moments. 5. Conclusions In this work, we have shown that DQMOM gives results equivalent to QMOM for the coagulation and sintering problems considered by Wright et al. (2001). In addition, we have shown that although all bivariate moments sets considered yield essentially the same results when N is large enough to achieve convergence, particular sets of moments are preferable for numerical reasons. Indeed, the real advantage of DQMOM over QMOM for the bivariate case is the ability to select arbitrary sets of independent moments without the need to find an algorithm to compute the weights and abscissas from the moments. The resulting well-conditioned set of ordinary differential equations can be rapidly solved using standard numerical methods, and the application of DQMOM to inhomogeneous cases is straightforward

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

1577

and does not directly solve for the moments. In contrast, spatial transport with QMOM is done in terms of the moments, and thus an algorithm is needed to compute the weights and abscissas from the moments at every time step. Finally, we have demonstrated that applying constraints to the source terms (such as forcing all weights to be constant in time) is generally not a good idea because it can lead to singularities caused by non-distinct abscissas. Appendix: Derivation of DQMOM equations The fundamental idea behind DQMOM is that we should choose the weights and abscissas such that as many moments computed from Eq. (9) as possible agree with moments computed from the transport equations (found from Eq. (3)). Note that there are a total of N weights and 2N abscissas and (equivalently) 3N unknown source terms in Eqs. (11)–(13). We will thus need to choose an equal number of independent moments to determine the source terms. The procedure for using these moments to find the source terms is described here. The space and time derivatives in Eq. (3) generate the corresponding terms in Eqs. (11)–(13). These are found by formally inserting Eq. (8), and differentiating (Fox, 2003; Marchisio & Fox, 2005). After some algebraic manipulations, this procedure yields jt n + u · jx n − jx · (Djx n) =

N


[(v − vn )(a − an ) + vn (1) (v − vn )(a − an ) + an (v − vn )(1) (a − an )]An

n=1

−

N


(1)

 (v − vn )(a

− an )Bn∗

−

n=1

−

N


N
n=1

(2) (v − vn )(a − an )wn Dvvn − 2

n=1

−

N


(v − vn )(1) (a − an )Cn∗ N


(1) (v − vn )(1) (a − an )wn Dvan

n=1

(v − vn )(2) (a − an )wn Daan ,

(29)

n=1

where An = jt wn + u · jx wn − jx · (Djx wn ),

(30)

Bn∗ = jt (wn vn ) + u · jx (wn vn ) − jx · (Djx wn vn ),

(31)

Cn∗ = jt (wn an ) + u · jx (wn an ) − jx · (Djx wn an ),

(32)

and

and the gradient-dependent functions appearing on the right-hand side are defined by Dvvn = D(jx vn ) · (jx vn ),

Dvan = D(jx vn ) · (jx an ),

and

Daan = D(jx an ) · (jx an ).

(33)

The derivatives of the delta function, denoted by (m) (x), have the property (Pope, 2000) 

+∞

−∞

(m) (x − s)g(x) dx = (−1)m g (m) (s),

where g (m) (x) is the mth derivative of g(x).

(34)

1578

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

We are now ready to apply the moment transform to Eq. (29). Formally, this yields  ∞ ∞ v k a l (jt n + u · jx n − jx · (Djx n)) dv da 0

0

=

N




l ∗ k l−1 ∗ (1 − k − l)vnk anl An + kv k−1 n an Bn + lv n an Cn − D(k, l),

(35)

n=1

where the gradient-dependent term is defined by D(k, l) =

N 


l−1 k l−2 k(k − 1)vnk−2 anl wn Dvvn + 2klv k−1 n an wn Dvan + l(l − 1)vn an wn Daan .

(36)

n=1

Comparing Eq. (35) to Eq. (10), we can see that formally jt Mkl + u · jx Mkl − jx · (Djx Mkl ) =

N 


l ∗ k l−1 ∗ (1 − k − l)vnk anl An + kv k−1 n an Bn + lv n an Cn − D(k, l).

(37)

n=1

For homogeneous systems, the terms involving spatial gradients will be null, and thus have no effect on the weights and abscissas. On the other hand, for inhomogeneous systems D(k, l) will be non-zero, and this term is needed to force the weights and abscissas to agree with the bivariate moments. One can note that in the absence of spatial diffusion (D = 0), D(k, l) will also be null. Finally, note that Eq. (37) is linear in the unknown source terms (An , Bn∗ , Cn∗ ). The next step is to consider the terms in Eq. (3) that correspond to transport in phase space. We begin by rewriting Eq. (3) as jt n + u · jx n − jx · (Djx n) = S,

(38)

where the phase-space transport term is defined by S ≡ −jv (Rv n) − ja · (Ra n) + Q.

(39)

We can then define the moment transform of the phase-space term by  ∞ ∞ S(k, l) ≡ v k a l S dv da. 0

(40)

0

As shown below, if S(k, l) is known, Eq. (37) forms a linear system that can be solved to find the unknown source terms. We can compute the phase-space moments using Eq. (39):  ∞ ∞ S(k, l) = − v k a l [jv (Rv n) + ja · (Ra n) − Q] dv da. (41) 0

0

As shown next, the integrals on the right-hand side can be approximated in terms of the weights and abscissas. Starting with the volume-growth term in Eq. (41), we can use integration by parts to find  ∞  ∞ k v jv (Rv n) dv = −k,0 Rv (0, a)n(0, a) − kv k−1 Rv (v, a)n dv, 0

(42)

0

where k,0 is the Kronecker delta. Using Eq. (8) to approximate n in the final integral, we find ∞ ∞

 0

0

k l

v a jv (Rv n) dv da = −k,0 Rv |v = 0nv (0) −

N


kwn vnk−1 anl Rv (vn , an ),

(43)

n=1

where Rv |v = 0 is the volume-growth rate of particles of zero volume and nv (0) is the corresponding number density function. Note that the first term on the right-hand side of Eq. (43) will be non-zero only for k = 0, and corresponds

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

1579

to the loss of particles due to evaporation (if Rv < 0) or to the gain of particles due to nucleation (if Rv > 0). In the absence of evaporation and nucleation, this term will be null. Turning next to the area-growth term in Eq. (41), we can use integration by parts to find  ∞  ∞ a l ja (Ra n) da = −l,0 Ra (v, 0)n(v, 0) − la l−1 Ra n da. (44) 0

0

Thus, the area-growth term becomes  0

∞ ∞ 0

v k a l ja (Ra n) dv da = −Ra |a = 0na (0) −

N


lw n vnk anl−1 Ra (vn , an ),

(45)

n=1

where Ra |a = 0 is the area-growth rate of particles of zero area and na (0) is the corresponding number density function. In general, this term will be null if amin in Eq. (7) is greater than zero. Turning now to the coagulation term, we will treat each of the two parts Q− and Q+ separately. The first part yields in a straightforward manner (Marchisio & Fox, 2005) ∞ ∞

 0

v k a l Q− dv da = −

0

N N



wn wq vnk anl
(vn , vq ).

(46)

n=1 q=1

The second part requires a change in the order of integration, and a change of variables (Marchisio & Fox, 2005; Wright et al., 2001). The final expression is  0

∞ ∞

v k a l Q+ dv da =

0

N N 1

wn wq (vn + vq )k (an + aq )l
(vn , vq ). 2

(47)

n=1 q=1

Note that the right-hand side of this expression is in closed form. Collecting together all of the terms, the phase-space moment term becomes S(k, l) = k,0 Rv |v = 0nv (0) + l,0 Ra |a = 0na (0) +

N 




l k l−1 kv k−1 n an wn Rv (vn , an ) + lv n an wn Ra (vn , an )

n=1

+

N N  1

wn wq (vn + vq )k (vn + vq )l − vnk anl − vqk aql
(vn , vq ). 2

(48)

n=1 q=1

Note that due to the form of the coagulation operator, the moments conserve volume (S(1, 0) = 0) and surface area (S(0, 1) = 0) when Rv and Ra are null. Thus, the weights and abscissas in the DQMOM representation will keep these conservation properties if the set of moments is chosen to include (k, l) = (1, 0) and (0, 1). Likewise, in the absence of coagulation and disappearance of particles due to evaporation, DQMOM will conserve number if (k, l) = (0, 0) is chosen. Hereinafter, we will neglect the source/sink terms appearing in Eq. (48) due to loss of particles of zero volume and/or area. The total particle number will then only decrease due to coagulation. From Eq. (10), we can now observe that the right-hand sides of Eqs. (37) and (48) must be equal. This equality provides a formal definition of An , Bn∗ and Cn∗ in terms of the weights and abscissas. Comparing the terms in Eqs. (37) and (48), we can note that the volume- and area-growth terms in the DQMOM representation can be solved for explicitly. This implies that any independent set of moments will yield exactly the same source terms for volume and area growth. Thus, the source terms can be written as Bn∗ = Bn + wn Rv (vn , an ),

(49)

Cn∗ = Cn + wn Ra (vn , an ),

(50)

1580

R.O. Fox / Aerosol Science 37 (2006) 1562 – 1580

where An , Bn and Cn are found from N 


l k l−1 (1 − k − l)vnk anl An + kv k−1 n an Bn + lv n an Cn = P (k, l) + D(k, l),

(51)

n=1

with the new term on the right-hand side depends only on coagulation: P (k, l) =

N N  1

wn wq (vn + vq )k (an + aq )l − vnk anl − vqk aql
(vn , vq ). 2

(52)

n=1 q=1

This is the form of the linear system used to find the source terms in the DQMOM transport equations appearing in the main text. References Diemer, R. B., & Ehrman, S. H. (2005). Pipeline agglomerator design as a model test case. Powder Technology, 156, 129–145. Diemer, R. B., & Olson, J. H. (2006). Bivariate moment methods for simultaneous coagulation, coalescence and breakup. Journal of Aerosol Science, 37, in press. Fan, R., Marchisio, D. L., & Fox, R. O. (2004). Application of the direct quadrature method of moments to polydisperse gas-solid fluidized beds. Powder Technology, 139, 7–20. Fox, R. O. (2003). Computational models for turbulent reacting flows. Cambridge: Cambridge University Press. Hindmarsh, A. C. (1983). ODEPACK, a systematized collection of ODE solvers. In Scientific computing. Amsterdam: North-Holland. Koch, W., & Friedlander, S. K. (1990). The effect of particle surface coalescence on the surface area of a coagulating aerosol. Journal of Colloid and Interface Science, 140, 419–427. Marchisio, D. L., & Fox, R. O. (2005). The direct quadrature method of moments for multivariate population balance equations. Journal of Aerosol Science, 36, 43–73. Marchisio, D. L., Pikturna, J., Fox, R. O., Vigil, R. D., & Barresi, A. A. (2003). Quadrature method of moments for population balances. AIChE Journal, 49, 1266–1276. Marchisio, D. L., Vigil, R. D., & Fox, R. O. (2003). Quadrature method of moments for aggregation-breakage processes. Journal of Colloid and Interface Science, 258, 322–334. McGraw, R. (1997). Description of aerosol dynamics by the quadrature method of moments. Aerosol Science and Technology, 27, 255–265. McGraw, R., & Wright, D. L. (2003). Chemically-resolved aerosol dynamics for internal mixtures by the quadrature method of moments. Journal of Aerosol Science, 34, 189–209. Piskunov, V. N., & Golubev, A. I. (1999). Treatment of size-dependent aerosol transport processes using quadrature based moment methods. Journal of Aerosol Science, 30, S451–S452. Piskunov, V. N., & Golubev, A. I. (2002). The generalized approximation method for modeling coagulation kinetics—Part 1: Justification and implementation of the method. Journal of Aerosol Science, 33, 51–62. Pope, S. B. (2000). Turbulent flows. Cambridge: Cambridge University Press. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (1992). Numerical recipes in Fortran 77. Cambridge: Cambridge University Press. Tandon, P., & Rosner, D. E. (1999). Monte Carlo simulation of particle aggregation and simultaneous restructuring. Journal of Colloid and Interface Science, 213, 273–286. Upadhyay, R. R., & Ezekoye, O. A. (2006). Treatment of size-dependent aerosol transport processes using quadrature based moment methods. Journal of Aerosol Science, 37, in press. Wang, L., & Fox, R. O. (2004). Comparison of micromixing models for CFD simulation of nanoparticle formation by reactive precipitation. AIChE Journal, 50, 2217–2232. Wright, D. L., McGraw, R., & Rosner, D. E. (2001). Bivariate extension of the quadrature method of moments for modeling simultaneous coagulation and sintering of particle populations. Journal of Colloid and Interface Science, 236, 242–251. Yoon, C., & McGraw, R. (2004a). Representation of generally mixed multivariate aerosols by the quadrature method of moments: I. Statistical foundation. Journal of Aerosol Science, 35, 561–576. Yoon, C., & McGraw, R. (2004b). Representation of generally mixed multivariate aerosols by the quadrature method of moments: II. Aerosol dynamics. Journal of Aerosol Science, 35, 577–598. Zucca, A., Marchisio, D. L., Barresi, A. A., & Fox, R. O. (2006). Implementation of the population balance equation in CFD codes for modelling soot formation in turbulent flames. Chemical Engineering Science, 61, 87–95.