Integrated Framework for the Numerical Solution of Multicomponent

4394

Ind. Eng. Chem. Res. 2004, 43, 4394-4404

Integrated Framework for the Numerical Solution of Multicomponent Population Balances. 2. The Split Composition Distribution Method Darren D. Obrigkeit, Timothy J. Resch, and Gregory J. McRae* Department of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139

The full multicomponent population dynamic equation is developed for systems with coagulation and solved using the split composition distribution method. A full discussion of the origins of this method and its advantages over conventional methods is followed by a complete and rigorous explanation of its numerical implementation. The Wiener expansion is introduced as a method for representing compositional variation in a number density distribution that is solved by orthogonal collocation on finite elements. Results are obtained by numerical integration using fifth-order Runge-Kutta and favorably compared to analytical results. Introduction Particle processing is a key step in numerous chemical production processes ranging from aerosols and inhaled drug delivery to protein characterization and structural genomics. Despite the widespread importance of these processes, our ability to accurately predict the performance of these systems is limited. First, the physics of these systems are complex, often involving multiple mechanisms acting over a wide range of particle sizes with several components undergoing simultaneous particle growth and coagulation. Second, current numerical methods are incapable of describing population balances necessary to model these systems with enough resolution to make meaningful comparisons with experimental data. These barriers have precluded thorough testing of the theoretical mechanisms used in simulations as well as the assimilation of data from experiments into process models. As a result, there is a fundamental lack of insight into the operation of these systems and a new generation of numerical models are required to meet the needs of designing, optimizing,1,2 and controlling3 these systems. Population balance problems are critical in multicomponent particulate systems. Typical phenomena that occur in these systems include coagulation, fragmentation, growth, condensation, and nucleation, among others. These phenomena are present in virtually every class of chemical processes including emulsions, copolymerization,4 hydrocarbons,4 and aerosols.5 Industrial applications range from nanostructured biomaterials to air pollution, pharmaceuticals, and mining. Despite the ubiquitous nature of particulate systems, accurate and computationally tractable models have not been developed because of the difficulty involved in framing the solution with a compact representation, especially for coagulation processes. This paper introduces a new representation for multicomponent particulate systems that is both accurate and computationally much faster than other methods and implements this representation to solve a two-component aerosol coagulation problem. Multicomponent population balances represent one of the most challenging sets of equations governing par* To whom correspondence should be addressed. Tel.: (617) 253-6564. Fax: (617) 258-0546. E-mail: [email protected].

ticle dynamics. As more components are added to these systems, the number of independent variables in the number density distribution n(m) increases. For a system with s components, the vector m contains s independent variables (m1, m2, ..., ms). The evolution of the distribution is described by the multicomponent population balance:6,7

|

∂n(m) ∂n(m) ) ∂t ∂t

+

coagulation

|

∂n(m) ∂t

|

∂n(m) ∂t

nucleation

+

+ fragmentation

|

∂n(m) ∂t

growth

+ ... (1)

As new components are included, new dimensions are added to the solution space for the population balance problem, and the solution time accordingly increases exponentially. In the case where p node points are used in each dimension, the solution time will scale as ps; this approach is generally prohibitive in any system that contains more than two or three components.8 Furthermore, when coagulation or fragmentation processes are active in a system, a number of integral relations must be incorporated into the general dynamic equation, transforming it from a partial differential equation (PDE) into a partial integrodifferential system. For multicomponent systems, the coagulation terms must be expressed as multidimensional integrals,9 which further contributes to the complexity of the problem and increases its computational intensity. As a result, new methods are needed to efficiently address these scaling issues to produce accurate, fast, and complete numerical solutions. Current efforts to this end have relied upon reductions in the multicomponent model to formulate a more compact representation and a more manageable set of equations for numerical solution. This work reviews multicomponent representations of population balance systems and describes current assumptions and methods for treating these problems. The split composition distribution is introduced as a new method that relaxes conventional constraints to permit solutions that more accurately reflect composition dynamics in a multicomponent system. The development, formulation, and numerical implementation of this method are described in detail, followed by a comparison with analytical results. Finally, the applicability of this

10.1021/ie0205492 CCC: $27.50 © 2004 American Chemical Society Published on Web 06/04/2004

Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4395 Table 1. Comparison of Solution Times as the Number of Components Increases no. of components

Kim and Seinfeld solution time (s)

split composition distribution solution time (s)

1 2 3 4

4 408 40540 59908

20 30 35 42

completely specify the particle characteristics. Referring to the third and forth rows of Table 2, we see that this can be accomplished either by defining the individual masses of each component in the particle

m ) (m1, m2, ..., ms)

(2)

or by defining the total particle mass and s - 1 mass fractions

(m, x) ) (m, x1, x2, ..., xs-1)

Figure 1. Batch crystallization initial and product distributions for a process optimized to produce the highest ratio of nucleated crystal mass/seed crystal mass possible. These results were reproduced using the method of characteristics and illustrate the presence of two distinct nodes in the product distribution due to seed and nucleated crystals.

method in conjunction with the node point placement methods introduced in part 1 of this series of papers is discussed. Representation of Population Balance Systems Multicomponent particle size distributions are described by multidimensional number density functions of the form n(O), where O represents a vector completely defining particle attributes. This number density is defined such that n(O) dφ1 dφ2 ... dφs represents the concentration of particles in the population with characteristics in the multidimensional hyperspace interval [φ1, φ1 + dφ1][φ2, φ2 + dφ2]...[φs, φs + dφs]. Depending on the system of interest, O can take on a number of different forms. Examples are given in Table 2 under the assumption that particle concentrations are given in terms of the number of particles per cubic meter. In a system such as the second row in Table 1, O takes on the defining population characteristic of particle mass, and the population distribution function is represented as in Figure 2. All of these results follow directly from early work regarding single-component coagulation and fragmentation pioneered in the 1960s.6,7 Multicomponent systems require that the independent variables

(3)

However, each of these systems have s components and therefore require ps equations if p is the number of node points along each independent axis. To produce a tractable problem, new representations must be introduced to reduce the equation set. One approach to reduce dimensionality involves using moments instead of discretized domains to describe the number density function.10,11 While this is convenient and effective at reducing solution times, moment methods produce very poor resolution results and cannot describe multimodal distributions that often occur as a result of coagulation and secondary nucleation in batch seeding and growth processes1 (see Figure 1). This represents a serious shortcoming because the exact shape of the particle size distribution is widely known to have significant effects on the processing efficiency2 and material properties. To reduce the dimensionality while preserving the resolution of the number density function, numerous techniques have been developed that simplify the variation of compositions while preserving the size resolution of the number density8,12-14 (see Figure 2). Several of these methods rely on discretization to address the difficulties associated with multicomponent particulate systems (see Table 3). The distribution splitting method12 proposes that multiple species be represented by a series of mixed- and pure-component distributions. For instance, a two-component system would be represented by four distributions: two distributions, each representing the number density of each pure component, another number density representing the mixed particles, and a final function describing the composition of the mixed particles. This method allows limited representation of variation in the composition at each size point in the overall population. For the two-component case, three different compositions are possible at any given mass point: pure component 1, pure component 2, and the composition given at that mass point in the mixed composition distribution. The solution time unfortunately scales as s × 2s-1 with the number of species s in the system. Another method that has been employed is the multicomponent sectionalization method.13 This method effectively discretizes the representation space along size and composition coordinates. While this method treats coagulation effects well, it does not conserve the number of particles in the system and poses significant difficulties in modeling all

Table 2. Multicomponent Number Density Distributions and Their Corresponding Cumulative Density Functions O v m

definition volume mass

n(O)

composition

no. of particles

1/(volume‚m3) n(v) dv

n(m)

1/(mass‚m3)

m ) {m1, m2, ..., ms} mass component n(m) {m, x} ) {m, x1, x2, ..., xs-1}

units

n(v)

1/(masss‚m3)

n(m,x) 1/(mass‚m3)

n(m) dm

cumulative particles

∫ n(v)dv N(m) ) ∫ n(m)dm N(m) ) ∫ ∫ ‚‚‚∫ n(m) dm dm ... dm N(m,x) ) ∫ ∫ ∫ ‚‚‚∫ n(m,x) dm dx dx

N(v) )

v

0

m

0

n(m) dm1 dm2 ... dms n(m,x) dm dx1 dx2 ... dxs

ms

0

m2

m1

0

xs-1

0

x2

0

0 x1

0

1

2

s

m

0

1

2

... dxs-1

4396 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004

Figure 2. Summary of multicomponent population balance representations and assumptions. Table 3. Comparison of Multicomponent Population Balance Representation Methods method full s-dimensional surface distribution splitting multicomponent sectionalization species mass distribution split composition distribution*

equation scaling ps ps2s-1 ps ps p(sn)a

no. of eqs, p ) 10, s ) 4, n ) 2

reference

10000 320 40 40 80

Kim and Seinfeld (1990) Toon et al. (1988) Gelbard and Seinfeld (1980) Pilinis (1990) This work

* a is related to the order of the polynomial chaos expansion and is 1 for most cases.

but the most elementary growth processes. The species mass distribution method14-16 addresses the number conservation caveats of the multicomponent sectionalization method by implementing a continuous analogue of the multicomponent sectionalization method. This can also be interpreted as reducing the dimensionality of multicomponent population balance problems via the internally mixed assumption, which states that particles of the same size all have the same composition17 and effectively reduces the governing equation set to ps equations representing each species in the system at p node points. The sacrifice for this reduction in complexity, of course, is the lack of ability to account for compositional variations at each size range of the population. In the past, the inability of measurement methods to simultaneously resolve size and composition variations within a population of particles has made the internally mixed assumption difficult to verify. Recent evidence, however, suggests that particles of the same size but very

different compositions coexist in atmospheric aerosols.5 Figure 2 summarizes the assumptions of these various methods and shows condensed versions of the coagulation contribution terms in the governing equations. An approach that can accommodate composition variations is proposed through the use of a Wiener expansion,18,19 also known as the polynomial chaos expansion, to compactly represent the composition at each size range. The representation of the multicomponent problem is now “split” into two separate distributions: (1) the number distribution as a function of the population member size and (2) the composition distribution at each size. This results in a method that scales linearly in s and produces solutions that are up to 105 times faster than conventional methods20,21 (see Table 1). For clarity, a complete and rigorous presentation of the assumptions and derivations for the split composition distribution method follows. Starting with the basic multicomponent population balance for coagulation

Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4397

∫0m ∫0m

∫0m β(m-c,c) n(m-c,t) n(c,t) ∞ ∞ ∞ dc1 ... dcs - n(m,t)∫0 ∫0 ‚‚‚∫0 β(m,c) n(c,t)

∂n(m,t) 1 ) ∂t 2

s

‚‚‚

s-1

1

dc1 ... dcs (4) the set of differential equations used in numerical implementations is developed. Results are presented and compared to a two-component analytical solution. Numerical Solution Methods for the Particle Size and Composition One method of representing a multicomponent mixture, known as the species mass distribution method, was introduced by Pilinis.14 This is a simplified approach that conveniently reduces the governing equation set with the assumption that all particles of the same size have the same composition, commonly known as the internally mixed assumption. While this assumption is not supported by experimental data,5 the resulting reduction in the equation set is necessary to produce a tractable solution time for conventional multicomponent population density representations. The species mass distribution is defined as qi(m,t), which represents the total mass of component i contained in particles size m, and is the product of the number density of particles size m with the mass of component i, mi, in each particle.

qi(m,t) ) mi(m,t) n(m,t)

(5)

Note that qi(m,t) dm is the total mass of component i in particles between sizes m and m + dm. Summing the individual species mass distributions and the individual species masses, mi, yields the total mass distribution and total mass, respectively:

q(m,t) ) m)

(6)

∂n(m) 1 ) ∂t 2

(7)

It follows that the contribution toward the mass of component i in particle size m from coagulation toward the mass distribution of particle size m is (see Figure 3)

∂qi(m,t) 1 ) ∂t 2

∫0mβ(c,m-c) (mi|c + mi|m-c)n(c,t) n(m-c,t) dc (8)

where β(c,m-c) is the coagulation kernel. Note that reducing the range of integration to (0, m/2) eliminates double counting of coagulation collisions and allows elimination of the factor 1/2.

∂qi(m,t) ) ∂t

∂qi(m,t) ) -mi|mn(m,t) ∂t

∫0∞β(m,c) n(c,t) dc

(10)

Together these two effects form the complete model for a system undergoing only coagulation:

∂qi(m,t) ) ∂t

∫0m/2β(c,m-c) (mi|c + mi|m-c)n(c,t) n(m-c,t) dc -

∫0∞β(m,c) n(c,t) dc

mi|mn(m,t)

(11)

qi(m,t) ) n(m,t) mi ) n(m,t) mxi(m,t)

(12)

Substituting this expression into eq 11 yields

Coagulation is the most difficult term to treat in population balance systems; reducing eq 4 produces the single-component coagulation expression:

∫0mβ(m-c,c) n(m-c,t) n(c,t) dc ∞ n(m,t)∫0 β(m,c) n(c,t) dc

The removal of component i in particle size m due to coagulation with other particles is

To implement this governing equation in a multicomponent model, we recall the relation

∑i qi(m,t) ∑i mi

Figure 3. Summary of multicomponent population balance representations and assumptions.

∫0m/2β(c,m-c) (mi|c + mi|m-c)n(c,t) n(m-c,t) dc (9)

∂[n(m,t) mxi(m,t)] ) ∂t

∫0m/2β(c,m-c) (mi|c + mi|m-c)n ∞ (c,t) n(m-c,t) dc - mi|mn(m,t)∫0 β(m,c) n(c,t) dc (13)

The split composition distribution method uses eq 13 in conjunction with the overall population balance, eq 9. The overall population balance determines the value of n(m), which is then used in conjunction with the solution of eq 13 to determine the value of the function xi(m,t). By choosing a representation for xi(m,t) that includes information about the variation of composition at each mass point m, the split composition distribution method solves the full multicomponent population balance, including composition variation among particles of the same size. This represents a significant advance over traditional methods such as the species mass distribution method,14 where particles of the same mass are restricted to having the same composition. This assumption implies that the product particle of mass m from any coagulation event would always have the same composition as the other particles of mass m. As shown in Figure 3, this assumption is clearly not reflective of coagulation processes and is therefore a key limitation of the species mass distribution representa-

4398 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004

Figure 4. Split composition distribution representation with the overall number density distribution and individual component distribution.

tion. When two particles of mass c and m - c coagulate, the mass of component i in the combined particle is mi|c + mi|m-c. The internally mixed assumption artificially forces this sum equal to mi|m, essentially contradicting the component mass balance. Split Composition Distribution Method To represent the composition distribution in a way that mirrors the true physics of population balance processes, an expansion developed by Wiener18 is employed. This expansion uses a set of orthogonal basis functionals of standard probability density to represent the distributions. In the case where independent Gaussian probability distributions of mean zero and standard deviation 1 [ξi(ω) ∼ N(0,1)] are used, Hermite polynomial functionals form the basis for the Wiener expansion:

xi ) xi,0 + xi,1ξ1(ω) + xi,2ξ2(ω) + xi,3[ξ1(ω)2 - 1] + xi,4ξ1(ω) ξ2(ω) + xi,5[ξ2(ω)2 - 1] + ... (14) where the xi,j’s are coefficients of the expansion. In this paper, Wiener expansions based on Gaussian probability density functions will be employed because of the convenient moment generating property, which enables automatic evaluation of the governing equations. The drawback of this basis function is a domain covering [-∞, ∞] instead of the natural [0, 1] domain for composition. Depending on the type of system represented and the required accuracy, different probability density functions may be used, varying number of terms included in the expansion. Recent applications of this technique in uncertainty analysis and representation of random variables discuss the selection of basis functions and convergence properties as more expansion terms are included.19,22 In the case of multicomponent population balance models, the Wiener expansion is used to describe the composition distribution of particles at each mass point in the population (see Figure 4). When a two-term Wiener expansion is used to represent the composition distribution of each component in the system

xi ) xi,0 + xi,1ξ1(ω)

(15)

then the composition distribution of each component i is Gaussian with mean xi,0 and standard deviation xi,1 (see Figure 5). Recalling that

mi|c ) cxi|c

(16)

we can then substitute the composition distribution

Figure 5. Two-term Wiener expansion representation for composition. This representation includes an overall number density function n(m) where the composition at each mass point m is described by a Wiener expansion. The first two terms in this expansion describe the mean and standard deviations of the composition of a particle size m.

expression from eq 15 into the governing equation for composition (eq 13) to produce

∂[n(m,t) m{xi,0(m,t) + xi,1(m,t) ξ1(ω)}] ) ∂t

∫0m/2β(c,m-c){c[xi,0(c,t) + xi,1(c,t) ξ1(ω)] +

(m - c)[xi,0(m-c,t) + xi,1(m-c,t) ξ1(ω)]}n(c,t) n(m-c,t) dc - m[xi,0(m,t) + xi,1(m,t) ξ1(ω)]n(m,t)

∫0∞β(m,c) n(c,t) dc

(17)

The complete set of governing equations for the multicomponent population balance are given by eqs 7 and 17. Equation 7 is a function of m and t and can therefore be numerically solved by discretizing along the m axis and performing numerical time step integration over the desired dynamic time range to produce the solution along the t axis. Equation 17 is a function of three variables, m, t, and ω. While m and t can be solved by discretization and time step integration, a solution along the ω axis must also be produced. This is accomplished as an error minimization problem where the error function R in eq 17 is reduced by orthogonalizing it with respect to the probability space ξi. The error function R represents the error as the difference between the right- and left-hand sides of eq 17:

R)

∂[n(m,t) m{xi,0(m,t) + xi,1(m,t) ξ1(ω)}] ∂t

∫0m/2β(c,m-c) {c[xi,0(c,t) + xi,1(c,t) ξ1(ω)] +

(m - c)[xi,0(m-c,t) + xi,1(m-c,t) ξ1(ω)]}n(c,t) n(m-c,t) dc + m[xi,0(m,t) + xi,1(m,t) ξ1(ω)]n(m,t)

∫0∞β(m,c) n(c,t) dc

(18)

The orthogonalization is completed by a set of normal equations that set the inner product of the error function R and each corresponding basis function of the Wiener expansion to zero:

Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4399

〈R, 1〉 ) 0 〈R, ξi(ω)〉 ) 0

(19)

Note that in this probability space the inner product is the integral of the inner product terms with the weighting function over the entire probability space:

〈R, f(m,t,ξi)〉 )

∫-∞∞Rf(m,t,ξi) pξ (ξi) dξi i

(20)

In this instance, the weighting function is a Gaussian probability density function pξ(ξi):

pξi(ξi) )

1 -ξ2i /2σ2 e σx2π

(21)

Conveniently, Gaussian probability density functions possess moment generating properties that simplify the evaluation of inner products:23

E[ξn] ) )

∫-∞∞pξ(ξ) ξn dξ ) ∫-∞∞

1 -ξ2/2σ2 n e ξ dξ x σ 2π

{

0 n ) 2k + 1 1 × 3 × ... (n - 1)σn n ) 2k

(22)

where k is any nonnegative integer. This allows easy reduction of the inner products given in eq 19 to produce the governing equations for composition:

∂[n(m,t) mxi,0(m,t)] m/2 ) 0 β(c,m-c) {cxi,0(c,t) + ∂t (m - c)xi,0(m-c,t)}n(c,t) n(m-c,t) dc -

∫

n(m,t) mxi,0(m,t)

∫0 β(m,c) n(c,t) dc

∂[n(m,t) mxi,1(m,t)] m/2 ) 0 β(c,m-c) {cxi,1(c,t) + ∂t (m - c)xi,1(m-c,t)}n(c,t) n(m-c,t) dc -

∫

n(m,t) mxi,1(m,t)

∫0 β(m,c) n(c,t) dc

(23)

Numerical Implementation

3

n(m)i ) Ri,1 + ηRi,2 + η Ri,3 + η Ri,4

(24)

Note that η is a local variable on each element and is equal to (m - mi-1)/(mi - mi-1)∈[0,1] within element i. Two conditions are imposed at the boundary of adjacent elements: (1) continuity and (2) continuity of the first

ηk

k

ηk

0

0

2

1 1 ≈ 0.789 + 2 x12

1

1 1 ≈ 0.211 2 x12

3

1

derivative. Continuity is an obvious requirement, and continuity of the first derivative imposes continuity of particle flux from particle growth between two elements, a necessary condition for any implementation that aspires to describe growth processes:

flux out of element i ) flux into element i + 1 -G

|

S

To numerically solve this system, the method of orthogonal collocation on finite elements is employed. This method converts the PDEs given by eqs 7 and 17 into ordinary differential equations (ODEs). By solving the resulting set of ODEs at a specially placed set of collocation points along the overall mass axis, the error in the solution is minimized and the optimal solution is obtained. Orthogonal collocation on finite elements uses a series of polynomial representation functions over local element domains illustrated in Figure 6. Typically, the polynomials for orthogonal collocation on finite elements are represented by spline curves. When a cubic-spline representation is chosen, n(m) is fully determined by four coefficients active on each element i: 2

Table 4. Collocation Point Placement for Cubic Splines over the Interval (0, 1) k

∞

∞

Figure 6. Placement of node points for orthogonal collocation on finite elements. The number density function n(m) is represented by a series of adjacent polynomial functions over elements defined by the boundaries m0, m1, m2, etc. At the boundaries of elements, adjacent polynomials are restricted to match the value and first derivative of n(m). The remaining degrees of freedom for each polynomial are determined by node points within each element.

|

∂n ∂n ) -G ∂m m∂m m+

(25)

where m- denotes the left hand side of the boundary located at the upper limit of element i and m+ denotes the right hand side of the boundary, the lower limit of element i + 1. Using the notation from eq 24, these two conditions are expressed as

Ri,1 + Ri,2 + Ri,3 + Ri,4 ) Ri+1,1 Ri,2 + 2Ri,3 + 3Ri,4 )

mi - mi-1 R mi+1 - mi i+1,2

(26)

The remaining two coefficients must be determined by two more conditions on each element. These conditions are determined by the value of n(m) at two collocation points ηk in each element:

4400 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 Table 5. Parameter Values Used in the Verification of the Two-Component Coagulation Case symbol

value

units

symbol

value

units

N0 b0

106 10-7

cm-3 cm3/s

m10 m20

10-13 10-15

g g

n(m)i,k ) Ri,1 + ηkRi,2 + ηk2Ri,3 + ηk3Ri,4

k ) 1, 2 (27)

Thus, for a system with N elements, there are 4N unknown coefficients, 2(N - 1) conditions imposed by eq 26 at the boundary between elements, and 2N more conditions imposed by eq 27. The remaining two conditions are determined by the value at each system boundary. The expressions in eq 26 are used to eliminate the third and fourth coefficients in eq 27 for every element but the last, producing a system of equations relating the coefficients to the values at the collocation points:

i ) 2, ..., N - 1; k ) 1, 2 i ) 1; k ) 0, 1, 2

n(m)i,k ) RN,1 + ηkRN,2 + ηk2RN,3 + ηk3RN,4 i ) N; k ) 1, 2, 3 (28) where the collocation points are placed at the roots of the second-order shifted Legendre polynomial on the interval η∈[0,1] ∼ m∈(mi-1, mi) corresponding to the interval covered by each element. This placement of collocation points minimizes the error in the solution and provides the optimal coefficients for the cubic-spline representation.24,25 Table 4 gives the location of collocation points for a cubic spline on each element. The points k ) 0 in the first element and k ) N in the last element are placed at the edges of the system boundary and serve as implementation points for the boundary conditions. Minimizing the error of the residual R for the differential equation at these collocation points produces the solution for the method of orthogonal collocation on finite elements. The residual is defined as the deviation between the left- and right-hand sides of the differential equation.

∂n(m,t) - f(m,t) dt

(29)

The cumulative error over the domain of interest (a, b) is represented as

error )

∫abR(m) dm

(30)

Gaussian quadrature is used to evaluate this integral exactly using a weighted sum of R(m) at specific mass points mi located at the roots of the corresponding orthogonal polynomial for the system (see Table 4):

error )

∫abR(m) dm ) ∑R(mi) wi i

∀ j ) 1, ..., P

(32)

where n(mj,t) is the number density at each of the P node points mj in the system. Note that the reduced system represented by eq 28 has 2N + 2 node points yi, k, and 2N + 2 unknown coefficients Ri,j. This system can be expressed in matrix form:

R ) A-1n(m)

(33)

where all position-dependent expressions are isolated in the A matrix, so the first and second derivatives with respect to m can be found by differentiating individual terms in the A matrix:

∂ ∂ ∂η ∂ n(m) ) ‚ n(m) ) A R ) BR ∂m ∂m ∂η ∂m

(

)

( )

∂2 ∂2 ∂η ∂η ∂2 ‚ ‚ n(m) ) n(m) ) A R ) CR ∂m ∂m ∂η2 ∂m2 ∂m2 (34) Resubstituting eq 33 yields two matrices, BA-1 and CA-1, which act as first and second derivative linear operators. These operators are convenient in evaluating derivatives such as those common in growth terms.

∂ n(m) ) BR ) BA-1n(m) ∂m ∂2 n(m) ) CR ) CA-1n(m) ∂m2

(35)

Equations 36-38 illustrate these matrix representations in full.

∂n(m,t) ) f(m,t) dt R(m) )

∂n(mj,t) ) f(m,t)|m)mj dt

n(m) ) AR

n(m)i,k ) (1 - 3ηk2 + 2ηk3)Ri,1 + (ηk - 2ηk2 + mi - mi-1 ηk3)Ri,2 + (3ηk2 - 2ηk3)Ri+1,1 + mi+1 - mi (ηk3 - ηk2)Ri+1,2

The error is minimized by forcing it to zero at each node point; this has the effect of forcing the error to be orthogonal to the cubic-spline solution representation, producing the most accurate solution for the given representation. This simply requires that R(mi) ) 0 for all mass points mi and is equivalent to requiring that the PDE given in eq 26 is satisfied at each node point. This produces a new system in the form

(31)

Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4401

Note that the terms 1/wi in the denominator are due to the chain rule and the fact that dη/dm ) 1/wi ) 1/(mi mi-1), where wi is the width of element i.

[

4402 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004

Orthogonal collocation on finite elements is used to discretize the m dimension for both the overall number density balance (eq 7) and the composition governing relationships given in eq 23. Because the composition distributions are described by spline functions of the composition mean and standard deviations, the composition distributions at each mass point in the distribution are easily reconstructed using standard Gaussian probability density functions.

∂g1(m1,m2) ∂g1(m1,m2) ∂m1 ∂m2 J(m1,m2) ) ∂g2(m1,m2) ∂g2(m1,m2) ∂m1 ∂m2

n(m,t) ) 4N0

s

1

∏ 2 i)1 m

J(m1,m2) )

(τ + 2)

i0

( ) [

exp -

∞

mi

1

∑ k)0

mi0

τ

s

mj

∏ s τ + 2 j)1 m

(k!)

j0

]

s

(41)

∏ i)1

e-mi /m20 mi0

(42)

For a binary system, this expression reduces to the first two terms in the product:

e-m1/m10 e-m2/m20 n(m1,m2) ) N0 m10 m20

(43)

However, a transformation is needed to translate the n(m1,m2) distribution into n(m,x1). The general formula for this transformation was previously derived in part 1 of this series of papers:

n(O) )

n(ψ) |det J(ψ)|

(44)

In this case, the vector O ) [m, x1] and ψ ) [m1, m2]. J(ψ) is given by the formula

Ji,j )

∂gi(ψ) ∂ψj

(45)

where gi(ψ) ) φi. For this application, g1(m1,m2) ) m and g2(m1,m2) ) x:

g1(m1,m2) ) m ) m1 + m2 g1(m1,m2) ) x ) The Jacobian J is given by

m1 m1 + m2

1 -m1

(m1 + m2)2 (m1 + m2)2

det J(m1,m2) ) -

m1 + m2 (m1 + m2)

2

)-

]

1 m

(48)

Taking the absolute value of the determinant produces the transformation

(40)

where β0 is the coagulation constant, N0 is the number of particles in the population, and t is time. The initial conditions for a binary system can be directly calculated from eq 40 for the case where τ ) 0:

n(m) ) N0

[

1 m2

k

where N0 is the initial number of particles in the population, mi0 is a parameter defining the initial distribution of component i, mj is the mass amount of component j, and τ is defined as

τ ) β0N0t

(47)

Evaluating the determinant of the Jacobian produces

Model Verification To verify the accuracy of the split composition distribution method, numerical results are compared against analytical results for a two-component exponential initial size distribution undergoing coagulation only. The analytical solution to this set of conditions is

]

(46)

n(m,x) )

n(m1,m2) |det J(m1,m2)|

) n(m1,m2)‚m

(49)

Using the n(m,x) distribution, it is possible to produce the overall number density distribution by integrating over all x values:

n(m) )

∫01n(m,x) dx

(50)

This analytical overall number density is compared to numerical results in Figure 7 for N0 ) 106, m10 ) 10-13, and m20 ) 10-15 at time points t ) 0, 20, 40, 100, 200, and 2000. Please refer to Table 5 for the full list of parameter values. Note that the numerical results closely match the analytical results over the full range of particle sizes and times, confirming the accuracy of this numerical method as well as the ability of this representation to describe the evolution of a particle size distribution. It is also possible to make a direct comparison of the average mass fraction of component 1 as a function of the particle size. The analytical mean mass fraction of component i is calculated as follows:

xji(m) )

total mass of component i in particle size m total mass of particle size m

∫01mxin(m,x) dx ∫01xin(m,x) dx ) ) n(m) ∫01mn(m,x) dx

(51)

Figures 8 and 9 compare the mean particle composition of the numerical and analytical solutions. The numerical results overall show good agreement with the analytical results, accurately reproducing trends across the entire particle size range. Conclusions The split composition distribution method offers a fast, accurate new approach for describing the variation of particle composition within a number density distribution. These methods scale linearly with the number of components in the system, producing solution times up to 3 orders of magnitude faster than conventional methods. The current implementation of this method chooses Gaussian probability density functions to de-

Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4403

Figure 7. Comparison of analytical and numerical results for a constant coagulation kernel operating on a two-component exponential initial particle size distribution. Over time, the particles represented in the leftmost particle size distribution coagulate to larger sizes and the distribution moves to the right. Note that the numerical results are almost in exact agreement with the analytical results over all ranges of particle size and time.

Figure 9. Close-up comparison of the numerical and analytical mean mass fraction x1. Note that the trends for this distribution are closely captured by the numerical solution, although they are not exactly reproduced because of the fact that only the first two moments of the composition distribution are used. The addition of more terms to the Weiner expansion could allow a more exact reproduction of the composition distribution.22

Acknowledgment This work was generously supported by the National Science Foundation, Bayer Corp., Aerodyne Research, Inc., and the National Center for Supercomputing Applications (NCSA). Nomenclature

Figure 8. Comparison of the numerical and analytical mean mass fractions x1 and x2.

scribe the composition variations within the population because of their unique moment generating properties. More terms can be added to increase the accuracy of the Wiener expansion.26 Although this requires more elaborate techniques to recover the composition distribution function, it essentially allows the user to customtailor the number of terms to their specific needs. For systems that include particle growth mechanisms, the Courant condition must be applied to ensure the stability of the solution. However, the same techniques for node point placement as those discussed in part 1 of this series of papers apply to the composition equations, allowing the user to optimize the node point placement and resolution over the particle size range, a technique that reduces the solution time by a factor of up to 107. In conjunction, the split composition distribution representation and growth mechanism node spacing methods enable the solution of many previously intractable problems. The flexibility of these techniques make them well-suited for a wide variety of critical applications ranging from pharmaceutical production to atmospheric pollution control.

E[x] ) expected value of the random variable x; calculated as the first moment of x with the probability density function for x as the weighting function J ) Jacobian operator, as defined in eq 45 n(m) ) number density of particles, where n(m) dm1 dm2 ... dms is the number of particles per unit volume in the multidimensional hyperspace interval [m1, m1 + dm1][m2, m2 + dm2]...[ms, ms + dms] mi|c ) mass of component i in a particle of total mass c pξi(ξi) ) probability density function for the random variable ξi(ω) qi(m) ) species mass balance for component i, representing the total mass of component i in particles of mass m Greek Symbols β(x,y) ) coagulation kernel (this function defines the rate at which particles of size x and y will coagulate to form a particle of size x + y) η ) local coordinate on a finite element, scaled linearly such that η ) 0 at the lower element boundary and η ) 1 at the upper element boundary ω ) “random act” that symbolizes the input for the “experiment” that produces the output of the corresponding random variable ξi(ω)27 ξi(ω) ) independent standard Gaussian random variable

Literature Cited (1) Chung, S. H.; Ma, D. L.; Braatz, R. D. Optimal Seeding in Batch Crystallization. Can. J. Chem. Eng. 1999, 77, 590. (2) Togkalidou, T.; Braatz, R. D.; Johnson, B. K.; Davidson, O.; Andrews, A. Experimental Design and Inferential Modeling in Pharmaceutical Crystallization. AIChE J. 2001, 47, 160. (3) Chiu, T.; Christofides, P. D. Nonlinear Control of Particulate Processes. AIChE J. 1999, 45, 1279. (4) Immanuel, C. D.; Cordiero, C. F.; Meadows, E. S.; Crowley, T. J.; Doyle, F. J. I. Modeling of Particle Size distribution in Emulsion Co-Polymerization: Comparison with Experimental Data and Parameteric Sensitivity Studies. Comput. Chem. Eng.

4404 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 (5) Noble, C. A.; Prather, K. A. Real-time Measurement of Correlated Size and Composition Profiles of Individual Atmospheric Aerosol Particles. Environ. Sci. Technol. 1996, 30, 2667. (6) Valentas, K. J.; Bilous, O.; Amundson, N. R. Analysis of Breakage in Dispersed Phase Systems. Ind. Eng. Chem. Fundam. 1966, 5, 271. (7) Valentas, K. J.; Amundson, N. R. Breakage and Coalescence in Dispersed Phase Systems. Ind. Eng. Chem. Fundam. 1966, 5, 533. (8) Kim, Y. P.; Seinfeld, J. H. Simulation of Multicomponent Aerosol Condensation by the Moving Sectional Method. J. Colloid Interface Sci. 1990, 135, 185. (9) Obrigkeit, D. D. Numerical Solution of Multicomponent Population Balance Systems with Applications to Particulate Processes. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, 2001. (10) Diemer, B. R.; Olson, J. H. A Moment Methodology for Coagulation and Breakage Problems: Part 1sAnalytical Solution of the Steady-State Population Balance. Chem. Eng. Sci. 2002, 57, 2183. (11) Diemer, B. R.; Olson, J. H. A Moment Methodology for Coagulation and Breakage Problems: Part 2sMoment Models and Distribution Reconstruction. Chem. Eng. Sci. 2002, 57, 2211. (12) Toon, O. B.; Turco, R. P.; Westphal, D.; Malone, R.; Liu, M. S. A Multidimensional Model for AerosolssDescription of Computational Analogs. J. Atmos. Sci. 1988, 45, 2123. (13) Gelbard, F.; Seinfeld, J. H. Simulation of Multicomponent Aerosol Dynamics. J. Colloid Interface Sci. 1980, 78, 485. (14) Pilinis, C. Derivation and Numerical-Solution of the Species Mass-Distribution Equations for Multicomponent Particulate Systems. Atmos. Environ. A: Gen. 1990, 24, 1923. (15) Katoshevski, D.; Seinfeld, J. H. Analytical Solution of the Multicomponent Aerosol General Dynamic Equationswithout Coagulation. Aerosol Sci. Technol. 1997, 27, 541. (16) Katoshevski, D.; Seinfeld, J. H. Analytical-Numerical Solution of the Multicomponent Aerosol General Dynamic Equationswith Coagulation. Aerosol Sci. Technol. 1997, 27, 550.

(17) Fernandez-Diaz, J. M.; Brana, M. A. R.; Garcia, B. A.; Muniz, C. G. P.; Nieto, P. J. G. Analytic Solution of the Aerosol Rigorous General Dynamic Equation without Coagulation in Multidimension. Aerosol Sci. Technol. 1999, 31, 3. (18) Wiener, N. The Homogeneous Chaos. Am. J. Math. 1938, 897. (19) Tatang, M. A. Direct Incorporation of Uncertainty in Chemical and Environmental Engineering Systems. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, 1995. (20) Kim, Y. P.; Seinfeld, J. H. Simulations of Multicomponent Aerosol Condensation by the Moving Sectional Method. J. Colloid Interface Sci. 1990, 135, 185. (21) Resch, T. J. A Framework for the Modeling of Suspended Multicomponent Particulate Systems with Applications to Atmospheric Aerosols. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, 1995. (22) Wang, C. Parametric Uncertainty Analysis for Complex Engineering Systems. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, 1999. (23) Papoulis, A.; Pillai, S. U. Probability, Random Variables, and Stochastic Processes; McGraw-Hill: Boston, MA, 2001. (24) Finlayson, B. A. Method of Weighted Residuals and Variational Principles; Academic Press: New York, 1972. (25) Villadsen, J.; Michelsen, M. L. Solution of Differential Equation Models by Polynomial Approximation; Prentice Hall: Englewood Cliffs, NJ, 1978. (26) Wang, M.; Liu, X. Y.; Sun, C.; Ming, N. B.; Bennema, P.; van Enckevort, W. J. P. Periodic Roughening Transitions in Diffusion-Limited Growth. Europhys. Lett. 1998, 41, 61. (27) Drake, A. W. Fundamentals of Applied Probability Theory; McGraw-Hill: New York, 1967.

Resubmitted for review January 14, 2004 Revised manuscript received January 14, 2004 Accepted March 1, 2004 IE0205492