Solution of Population Balance Equations for Emulsion Polymerization

Ind. Eng. Chem. Res. 2007, 46, 643-654

643

Solution of Population Balance Equations for Emulsion Polymerization: Zero-One and Zero-One-Two Systems Hugo M. Vale and Timothy F. McKenna* CNRS-LCPP/ESCPE-Lyon, BP 2077, 43 BlVd du 11 NoVembre 1918, 69616 Villeurbanne Cedex, France

A numerical method was developed to solve the sets of population balance equations (PBEs) that govern the evolution of the particle size distribution (PSD) in emulsion polymerization processes and, more specifically, in zero-one and zero-one-two systems. The technique is based on the finite-volume method and uses high-resolution schemes and a generalized fixed pivot technique to discretize the growth and aggregation terms, respectively. Implicit-explicit methods are used to integrate the semidiscrete equations in time. The performance of the method was evaluated by comparing the computed distributions and/or their moments to analytical solutions determined for selected cases. In all four cases analyzed, good qualitative and quantitative agreement was found between numerical and analytical results. The present algorithm made possible to simulate, for the first time, the evolution of PSD for an actual emulsion polymerization system using the zero-onetwo model. 1. Introduction The modeling of particle size distribution (PSD) in emulsion polymerization is a subject of particular interest both from a fundamental (mechanistic studies) and an applied (reactor modeling) perspective. A thorough discussion on this topic can be found in a review1 recently published by the present authors. The dynamic evolution of the latex PSD is mathematically described by a population balance equation (PBE), or a system of PBEs, that accounts for the various phenomena capable of affecting particle size and number; namely, growth by polymerization, particle coagulation, particle nucleation, and radical kinetics. The precise form of the population balance model depends on the assumptions made in treating radical kinetics, which in turn are usually dictated by the characteristics of the system and the objectives pursued. The simplest approach, and also the most commonly used in polymer reaction engineering, is the pseudo-bulk (PB) model, which neglects compartmentalization between particles of the same size.1 Another option, more complex and particularly used in mechanistic studies, is the zero-one (0-1) model, which assumes that particles contain at most one radical.1,2 An even more sophisticated possibility is the zero-one-two3-5 (0-1-2) model, which distinguishes between particles having zero, one, or two radicals. The resolution of PBEs is, in general, a difficult and problemspecific task. For a general survey of techniques, the reader is referred to refs 6-10. Numerical methods applied up to now to solve emulsion polymerization PBEs have been critically reviewed by Vale and McKenna.1 As pointed out by these authors, a particularly efficient and easy-to-implement approach consists of using the finite-volume (FV) method combined with a semidiscrete high-resolution11 scheme to approximate the hyperbolic (growth) term and with a suitable discretization formula for the coagulation terms, such as the well-reputed fixed pivot technique.12 Lim et al.13 were apparently the first to employ such a combination of techniques to the solution of PBEs (for crystallization), though they actually reported erroneous expressions for the aggregation terms. The merits of the approach just described are multiple: * To whom correspondence should be addressed. Tel.: +33-4-72431775. Fax: +33-4-724-31768. E-mail: [email protected].

1. The FV formulation guarantees that the number of particles is conserved.11 2. The numerical difficulties related to the advective term (dissipation and/or dispersion) are subdued by the use of a highresolution scheme. Thus, unlike other methods,9,14 there is no need to add an artificial viscosity term to eliminate spurious oscillations. 3. Singular sources such as particle nucleation can be easily treated, in contrast with other methods10,15 that require special stratagems. 4. The internal consistency of the fixed pivot technique ensures that the coagulation terms are approximated in a way that preserves particle mass and number. 5. It is applicable to arbitrary nonuniform grids. The aim of the present work is to extend this approach, which so far has only been applied1 to the PB model (single PBE), to solve the intricate coupled systems of PBEs that arise in 0-1 and 0-1-2 models. Thus far, the only solutions obtained for the full 0-1 model (i.e., including coagulation) have been obtained with classical finite-difference (FD) methods,16 which are not appropriate for such problems.1 With respect to the full 0-1-2 model, to the best of our knowledge, no solutions have been obtained. The organization of the paper is as follows. In Section 2, we present the PBEs for the 0-1 and 0-1-2 models. Then, in Section 3, the numerical method is described in detail. The performance of the method is evaluated in Section 4, by comparing numerical results against analytical solutions determined for selected cases. Concluding remarks are given in Section 5. 2. Population Balance Equations There are two alternative formulations of the 0-1 model: one that just distinguishes between particles having zero or one radicals and another that further distinguishes between monomeric or polymeric radicals.1,2 To differentiate them, hereafter these will be symbolized 0-1-S and 0-1-MP, respectively. With regard to the 0-1-2 model, the formulation analyzed in this work does not account for radical chain length, as is the case in virtually all reactor models; more detailed approaches5 are still out of reach. For the sake of completeness and for

10.1021/ie060928l CCC: $37.00 © 2007 American Chemical Society Published on Web 12/22/2006

644

Ind. Eng. Chem. Res., Vol. 46, No. 2, 2007

Figure 1. Discretization of the unswollen radius domain.

The discretization of the size coordinate is depicted in Figure 1. A cell-centered mesh is built by dividing the truncated radius domain [rnuc,rmax] into M contiguous cells Ci ) [ri-1/2,ri+1/2], 1 e i e M, of size ∆ri ) (ri+1/2 - ri-1/2) and center ri ) 1/2(ri-1/2 + ri+1/2). The average value of f(r,t) over Ci at time t is denoted

Table 1. Emulsion Polymerization Models and Corresponding Distributions model

distributions f)f f ) f ) f )

PB 0-1-S 0-1-MP 0-1-2

[f0 f1 ]T [f0 f1M f1P ]T [f0 f1 f2 ]T

∴

f1 ) f1M + f1P

comparison, the equations for the PB model are also included in this presentation. Since we will be dealing with multiple PBEs, it is convenient to define a vector distribution f(V,t) that regroups the individual number density functions (Table 1). The overall distribution is given by the sum of the individual distributions, f(V,t) ) ∑mfm(V,t). The reader is referred to the Nomenclature section for the meaning of the variables. For an ideal stirred tank reactor, all population balance models considered here can be written in the standard form of a system of linear hyperbolic balance laws

∂ ∂ f(V,t) ) - [v3 (V,y)f(V,t)] + s(f,V,y) ∂t ∂V

(1)

where v3 (V,y) is a row vector of volume growth rates, s(f,V,y) is the source function, and y(t) is the continuous-phase vector. The source function is given by

s(f,V,y) ) R(V,y)f(V,t) + a(f,V,y) +

(

)

1 in f (V,t) θin

V˙ w 1 + f(V,t) + nδ(V - Vnuc)Rnuc(y) (2) out Vw θ

where the various terms account, respectively, for kinetic coupling between the different types of particles,17 particle aggregation, particle inflow, particle outflow plus volume change, and particle nucleation. The definitions of v3 , R, a, and n for the various models are summarized in Table 2. Note that, by summing up along the columns, the 0-1 and 0-1-2 equations reduce to the PB equation.

The method proposed here is of the semidiscrete type. The equations are first discretized only in space (here, “space” refers to the internal particle size coordinate), resulting in a system of ordinary differential equations (ODEs) in time, which is subsequently integrated using an appropriate ODE method. 3.1. Spatial Discretization. In the previous section, the PBEs were stated in terms of the unswollen volume. This choice simplifies the writing of the aggregation terms and facilitates the derivation of analytical solutions. On the other hand, numerical solutions are preferably obtained using the unswollen radius as the internal coordinate. In terms of r, eq 1 takes the form

with r3 (r,y) ) v3 (V,y)/(4πr2) and fr(r,t) ) 4πr2fv(V,t).

∫rr

1 ∆ri

i+1/2

(3)

f(r,t) dr

(4)

i-1/2

Integrating eq 3 over Ci and dividing by ∆ri, we obtain the following system of semidiscrete equations (1 e i e M)

1 d [h(fi+1/2,ri+1/2,y) - h(fi-1/2,ri-1/2,y)] + jsi(f,y) hf (t) ) dt i ∆ri (5) where h(f,r,y) ) r3 (r,y)f(r,t) is the (variable coefficient) flux function and jsi is the average value of the source function over Ci. Note that eq 5 is, up to this point, an exact relation. To construct a numerical method, we must define suitable approximations to h(fi+1/2,ri+1/2,y) and jsi. Appropriate numerical fluxes can be computed by a variety of high-resolution schemes.11,18 In this work, we employ the widely used finite-volume WENO (weighted essentially nonoscillatory) schemes,19,20 but obviously other methods may be utilized if desired. More specifically, we employ a fifth-order formula. Detailed instructions for the implementation of WENO algorithms can be found in the notes of Shu20 and, thus, will not be repeated in this text. Here, it suffices to say that for a system such as the one at hand, where r3 g 0, the numerical flux is simply given by hi+1/2 ) r3 (ri+1/2,y)fi+1/2 . The value fi+1/2 is the left-side approximation to the pointwise value of f(r,t) at the (i + 1/2)th cell boundary (see Figure 1) and is obtained from the known cell averages hfi(t) by means of the WENO reconstruction procedure. Let us now address the source term. Integrating the expression for s(f,r,y) over Ci and dividing by ∆ri, we obtain

1 in hf i (t) θin V˙ w δi,1 1 + hf i(t) + n Rnuc(y) (6) out V ∆r θ w i

jsi(f,y) ) [R(r,y)f(r,t)]i + aji(f,y) +

3. Numerical Method

∂ ∂ f(r,t) ) - [r3 (r,y)f(r,t)] + s(f,r,y) ∂t ∂r

hf i(t) )

(

)

where we have made use of the fact that nucleated particles enter the first cell. Approximations must be formulated for the first two terms on the RHS of eq 6. For the kinetic coupling term, we have opted for the simplest solution, which consists of applying the second-order approximation

[R(r,y)f(r,t)]i ) R(ri,y)fhi(t) + O(∆ri2)

(7)

Higher-order approximations could be obtained by means of an appropriate quadrature rule, but at the expense of more CPU time (to reconstruct f at the corresponding integration nodes). Besides, as shown in Section 3.3, this formulation is convenient for integration over time.

Ind. Eng. Chem. Res., Vol. 46, No. 2, 2007 645 Table 2. Definition of the Terms Appearing in Equation 1 for the Various Types of Modelsa

a

For simplicity, dependence on t, V, and y is omitted when possible.

It remains to develop a proper approximation for the aggregation term, which is the subject of the following section. 3.2. Discretization of the Aggregation Integrals. The discretization of the aggregation integrals appearing in PBEs is a delicate task. Indeed, special measures must be taken so that the method does not produce physically incoherent results (e.g., nonconservation of mass upon aggregation). The solution we propose in this work is to generalize the well-reputed fixed pivot technique,12 originally derived for a single (scalar) PBE, to a system of PBEs. This method guarantees preservation of particle number and mass and is compatible with the finite-volume formulation.

In terms of volume distributions, the aggregation vector, a, is composed of birth and death terms whose generic formulas are (cf. Table 2)

Bv ) (1 - 1/2δmq)

∫V

V-Vnuc nuc

β(V - V′,V′)fm(V - V′)fq(V′) dV′

Dv ) -fm(V)

∫V∞

β(V,V′)f (V′) dV′

(8)

(9)

nuc

where the superscript specifies the internal coordinate. By applying the extended fixed pivot technique to eqs 8 and 9, we will obtain approximations for the cell averages B h iv and D h iv. These can then be converted into their radial equivalents by means of the trivial relations B h ir∆ri ) B h iv∆Vi and D h ir∆ri ) D h iv∆Vi.

646


For the reader less acquainted with the method of Kumar and Ramkrishna,12 we recapitulate here its basic principles: (i) the particles in each cell are assumed to be concentrated at a representative size, named the pivot; (ii) if the size of a particle formed by aggregation does not fall on any pivot, the new particle is assigned to the two nearby pivots; (iii) the fraction attributed to each pivot is determined by imposing the conservation of two properties. In this work, the pivots are assumed to be the center of the cells, that is, the points ri (or in terms of volume, Vi ) (4/3)πri3). Accordingly, the ith cell will receive a contribution for every particle born in the region Ri ) [Vi-1,Vi+1]. For the preservation of number and mass, this fraction reads12

{

Vˆ - Vi-1 , Vi-1 e Vˆ eVi Vi - Vi-1 ηi(Vˆ ) ) V - Vˆ i+1 , V e Vˆ eVi+1 Vi+1 - Vi i

(10)

where Vˆ is the unswollen volume of the new particle. Since the particles in Ci are assumed to be concentrated at the pivotal point, we express the number density functions as M

fm(V,t) )

N mi (t)δ(V - Vi) ∑ i)1

(11a)

M

f(V,t) )

Ni(t)δ(V - Vi) ∑ i)1

(11b)

where Nmi and Ni ) ∑mN mi are, respectively, the number of particles of type m and the total number of particles in cell i. These particle numbers are accessible from the average cell values of the corresponding number density functions (cf. eq 4). We can now move on to the discretization of eqs 8 and 9, starting by the birth term. Recalling that Ci receives a fraction ηi for every particle formed in the region Ri, we can integrate eq 8 to obtain

∫C Bv dV ) (1 - 1/2δmq)∫R dV ηi(V)∫VV-V i

i

nuc

Replacing fm and fq according to eq 11a, dividing by ∆ri, and after some manipulations to process the Dirac deltas, we get as final expressions (see Appendix A)

ηi(Vˆ )β(rj,rk,y)N mj (t)N qk (t), ∑ ∆r j,k∈Q i

Qi ) {{j,k}: 1 e k,j e i; Vi-1 < Vˆ < Vi+1} B h ri

1

)

∑

∆ri j,k∈Qi

(1 -

(13)

where Vˆ ) Vj + Vk. An analogous treatment is used for the death term. We first integrate eq 9 over Ci to obtain

∫C

i

∫

i

∫V

∞

β(V,V′)f(V′) dV′

∑β(ri,rk,y)Nk(t)

(16)

k)1

)

(17)

where LE and LI denote, respectively, the components of u˘ (t) to be treated explicitly and implicitly. This method was selected for its well-known favorable stability properties. For the first step, which cannot be done with eq 17, we use the forwardbackward Euler scheme (IMEX-BDF1)18

(15)

nuc

and then replace fm and f according to eq 11, divide by ∆ri, and

(18)

To apply the above schemes to eq 5, we start by defining a vector that regroups all terms to be treated explicitly

ei ) -

1/2δjk)ηi(Vˆ )β(rj,rk,y)Nmj (t)Nmk (t)

dV fm(V) C

∆ri

un+1 ) un + ∆t(LEn + LIn+1)

Qi ) {{j,k}: 1 e k e i; k e j e i; Vi-1 < Vˆ < Vi+1} (14)

Dv dV ) -

M

Nmi (t)

Equations 13, 14, and 16 constitute the set of relations necessary to compute the cell average aji appearing in eq 6. The possible combinations of cells j and k leading to the formation of a particle in the region Ri (eqs 13 and 14) are computed just once, at the moment of grid generation. Likewise, the values of ηi(Vj + Vk), which are time independent, need only to be calculated once. Computationally, efficient storage of these variables can be achieved through the use of pointers. 3.3. Time Integration. The ODE system represented by eq 5 must be integrated in time, together with the differential and/or algebraic equations for y(t). In general, a standard ODE solver cannot be used to integrate eq 5 (this was confirmed), since it involves terms of different types. The kinetic coupling term, Rf, is very stiff and therefore demands an implicit method. The advective (growth) and aggregation terms, on the other hand, are preferably treated via explicit methods as they are nonlinear and highly coupled over the species and over space. With respect to the advective term, it is worth stressing that even though implicit schemes for advection allow larger Courant numbers as far as stability is concerned (which would be very attractive), with regard to positivity and TVD they are only marginally superior to explicit methods.18 Hence, there is no advantage in treating the advective term implicitly. The way to deal with the distinct nature of the terms is to treat them separately, each with an optimal method. Specifically, we would like to treat the kinetic coupling term with an implicit scheme and the remaining terms with an explicit scheme. Two approaches are widely used to do this: operator splitting11,18 and IMEX18 (implicit-explicit) methods. In this work, we use the latter approach, given its advantages18 over the first. Among the various IMEX methods available in the literature, we have adopted the second-order IMEX-BDF2 scheme18,21

m*q

i

1

(

nuc

1

D h ri ) -

4 1 2 2 4 un+1 ) un - un-1 + ∆t LEn - LEn-1 + LIn+1 3 3 3 3 3

β(V - V′,V′)

fm(V - V′)fq(V′) dV′ (12)

B h ri )

process the Dirac deltas. As shown in Appendix A, the final result is

(

)

V˙ w 1 1 1 (hi+1/2 - hi-1/2) + aji + inhf iin - out + hf + ∆ri Vw i θ θ δi,1 n Rnuc (19) ∆ri

The time discretization of eq 5 on the basis of the IMEX-BDF2 scheme then reads

4 4 1 2 hf in+1 ) hf in - hf in-1 + ∆t ein - ein-1 + 3 3 3 3 2 ∆tR(ri,yn+1)fhin+1 (20) 3

(

)

Ind. Eng. Chem. Res., Vol. 46, No. 2, 2007 647

Because the implicit component is actually linear in f, thisexpression can be easily rearranged to obtain a relationship explicit in hfin+1

2 hf in+1 ) I - ∆tR(ri,yn+1) 3

[

-1

4 n 1 n-1 4 hf - hf i + ∆t ein 3 i 3 3 2 n-1 (21) e 3 i

][

(

)]

Table 3. Dimensionless PBEs

[] [][

/ -R ∂ f0 ∂ 0 )+ ∂τ f /1 ∂V˜ f /1 R

[

(22)

The inverse matrices appearing in the last two equations are calculated analytically as functions of the coefficients of R, so as to avoid repeated numerical inversions. In addition, note that the schemes require R(ri,yn+1), which is generally not available at t ) tn(unless R is time independent). An iterative procedure is therefore necessary. We start with R(ri,yn) to obtain a first estimate of hfin+1, which is then used to compute an approximation of yn+1 (by solving the corresponding differential algebraic equations (DAEs)). A new estimate of R(ri,yn+1) is determined, and the cycle continues until convergence is attained. Typically, the number of iterations required is small and the CPU time involved is unimportant with respect to that of evaluating ein. The time integration is performed with a constant step size, chosen small enough such that the stability of the method (dictated by the hyperbolic term) is ensured and that numerical inaccuracies are mainly due to the spatial discretization. Variable step-size versions of IMEX multistep methods also exist,22 which in principle permit the implementation of automatic stepsize control. Nevertheless, the computation of suitable error estimates for such methods is still the object of current research. 4. Numerical Results In this section, the proposed numerical algorithm is validated against certain analytical solutions. In addition, results gained from the simulation of an actual emulsion polymerization system are presented to illustrate the capabilities of the method in resolving practical problems. Few analytical solutions have been reported in the literature for the systems addressed in this work. To the best of our knowledge, the only solutions refer to the 0-1-S model and are restricted to the particular cases of seed growth,17,23 and simultaneous nucleation and growth.17,24 To minimize this gap, additional solutions were derived. Although the set of analytical results gathered here is still limited, it should nevertheless be sufficient to assess the capacity of the numerical method to deal with the most problematic terms of eq 3: growth, kinetic coupling, nucleation, and aggregation. In deriving new solutions, we concentrated our efforts in the 0-1-S and 0-1-2 models, because the number of parameters is smaller and thus analytical solutions are expected to be less cumbersome. The following assumptions are made: (i) closed system (θ f ∞), (ii) constant reaction volume (V˙ w ) 0), (iii) constant coefficients, (iv) size-independent aggregation kernel (β ) β0), and (v) Vnuc ) 0. Additionally, the following dimensionless quantities are defined: f /m(V˜ ) ) fm(V)V0/N0, V˜ ) V/V0, τ ) Kt/V0, R ) FV0/K, κ ) kdes/F, γ ) c/F, Λ ) β0N0V0/K, and Ω ) RnucV0/(KN0). Because numerical solutions are obtained in terms of r, we further define f /m(r˜) ) fm(r)r0/N0 and r˜ ) r/r0, from which it follows that f /m(r˜) ) 3r˜2f /m(V˜ ). Under the aforesaid assumptions, the PBEs for the 0-1-S and 0-1-2

R(1 + κ) -R(1 + κ)

][ ] f /0

f /1

+

∫ f (V˜ - V˜ ′)f (V˜ ′) dV˜ ′ + 1/2∫ f (V˜ - V˜ ′)f (V˜ ′) dV˜ ′- f (V˜ )∫ f*(V˜ ′) dV˜ ′ ∫ f (V˜ - V˜ ′)f (V˜ ′) dV˜ ′ - f (V˜ )∫ f*(V˜ ′) dV˜ ′ V˜

1/2

0

Λ

Proceeding in an analogous manner with the IMEX-BDF1 scheme, we get

hf in+1 ) [I - ∆tR(ri,yn+1)]-1(fhin + ∆tein)

0-1-S Model

/ 0

V˜

/ 0

V˜

0

/ 1

0

/ 0

/ 1

/ 1

∞

/ 0

0

∞

/ 1

0

[] [ ][

[

/ 0 / 1 / 2

0 -R ∂ f/ )+ R 1 ∂V˜ 0 2f /2

0-1-2 Model Rκ -R(1 + κ) R

2Rγ R(1 + 2κ) -R(1 + 2κ + 2γ)

0 δ(V˜ )

][ ] f /0

f /1 + f /2

∫ f (V˜ - V˜ ′)f (V˜ ′) dV˜ ′ + 1/2∫ f (V˜ - V˜ ′)f (V˜ ′) dV˜ ′ - f (V˜ )∫ f*(V˜ ′) dV˜ ′ ∫ f (V˜ - V˜ ′)f (V˜ ′) dV˜ ′ +∫ f (V˜ - V˜ ′)f (V˜ ′) dV˜ ′ - f (V˜ )∫ f*(V˜ ′) dV˜ ′ ∫ f (V˜ - V˜ ′)f (V˜ ′) dV˜ ′ + 1/2∫ f (V˜ - V˜ ′)f (V˜ ′) dV˜ ′ - f (V˜ )∫ f*(V˜ ′) dV˜ ′ V˜

1/2

Λ

+

[ ]

Ω

f ∂ f ∂τ f

]

0

V˜

0 V˜

0

/ 0

V˜

/ 0

/ 0

/ 1

/ 0

/ 2

0

V˜

0

/ 1

/ 2

/ 2

V˜

0

/ 2

/ 1

∞

/ 0

∞

/ 1

/ 1

0

0

/ 2

∞

0

]

[ ]

0 + Ω δ(V˜ ) 0

Table 4. Moment Equations for the PBEs of Table 3 0-1-S Model λ˙ 0 ) -Rλ0 + R(1 + κ)λ1 + Λ[(λ02 + λ12)/2 - λ0λ] λ˙ 1 ) Rλ0 - R(1 + κ)λ1 + Λ(λ0λ1 - λ1λ) + Ω µ˘ 0 ) -Rµ0 + R(1 + κ)µ1 + Λ(λ0µ0 + λ1µ1 - λµ0) µ˘ 1 ) Rµ0 - R(1 + κ)µ1 + Λ(λ1µ0 + λ0µ1 - λµ1) + λ1 0-1-2 Model λ˙ 0 ) -Rλ0 + Rκλ1 + 2Rγλ2 + Λ[(λ02 + λ22)/2 - λ0λ] λ˙ 1 ) Rλ0 - R(1 + κ)λ1 + R(1 + 2κ)λ2 + Λ(λ0λ1 + λ1λ2 - λ1λ) + Ω λ˙ 2 ) Rλ1 - R(1 + 2κ + 2γ)λ2 + Λ(λ0λ2 + 1/2λ12 - λ2λ) µ˘ 0 ) -Rµ0 + Rκµ1 + 2Rγµ2 + Λ(λ0µ0 + λ2µ2 - µ0λ) µ˘ 1 ) Rµ0 - R(1 + κ)µ1 + R(1 + 2κ)µ2 + Λ(λ1µ0 + λ0µ1 + λ2µ1 + λ1µ2 - µ1λ) + λ1 µ˘ 2 ) Rµ1 - R(1 + 2κ + 2γ)µ2 + Λ(λ2µ0 + λ0µ2 + λ1µ1 - µ2λ) + 2λ2

models (cf. Table 2) can be rewritten in dimensionless form as shown in Table 3. The equations that describe the evolution of the number of particles having m radicals

λm(τ) )

∫0∞f /m(V˜ ,τ) dV˜ ) ∫0∞f /m(r˜,τ) dr˜

(23)

and the volume of particles having m radicals

µm(τ) )

∫0∞ V˜ f /m(V˜ ,τ) dV˜ ) ∫0∞r˜3f /m(r˜,τ) dr˜

(24)

may be obtained by standard procedures6 from the equations of Table 3 and are summarized in Table 4. At any point in time, the total number of particles and the total volume of particles are given, respectively, by λ ) ∑mλm and µ ) ∑mµm. Despite the simplifications made, finding solutions to the PBEs of Table 3 (via Laplace transforms) can only be done in special cases. Table 5 indicates the cases considered and the types of solutions that are possible to obtain. Often, the Laplace transforms of the solutions are too complicated for analytical inversion. These transforms can, nevertheless, be easily and accurately inverted by a suitable algorithm,25 and therefore, for

648


Table 5. Summary of Cases Considered and Solutions Obtaineda 0-1-S case 1 2 3 4

f seed growth Ω)Λ)0 nucleation and growth Λ)0 aggregation Λ f ∞, τ ) N0β0t growth, nucleation, and aggregation

/ m

λm

0-1-2 µm

L

A

A

L

A

A

A

A

A

A/N

N

f

/ m

λm

µm

A

N

L

A

N

A/L

A

A

A/N

N

a A: analytical solution. L: Laplace transform of the solution. N: numerical solution.

Table 6. Results for Case 1 τ 10

50

100

0-1-S: R ) 1, κ ) 1, M ) 100 0.05 0.11 L2 error f /m (%) error λm (%) 0.00 0.00 error µm (%) 0.01 0.01

0.17 0.00 0.01

L2 error f /m (%) error λm (%) error µm (%)

0-1-S: R ) 10, κ ) 1, M ) 100 0.56 2.52 0.00 0.00 -0.16 -0.67

4.13 0.00 -0.81


0-1-S: R ) 10, κ ) 1, M ) 200 0.05 0.28 0.00 0.00 0.00 -0.04

0.62 0.00 -0.08

0-1-2: R ) 1, κ ) 2.11, γ ) 2, M ) 100 error λm (%) 0.00 0.00 error µm (%) 0.01 0.01

0.00 0.01

0-1-2: R ) 10, κ ) 2.11, γ ) 2, M ) 200 error λm (%) 0.00 0.00 error µm (%) 0.00 -0.03

0.00 -0.06

practical purposes, solutions obtained by this method may be taken as analytical. Note also that moment solutions can always be obtained, either analytically or numerically, which guarantees us at least an element of comparison for every case. Appendix B regroups all analytical solutions derived. 4.1. Case 1. We first consider the case in which only growth by polymerization takes place. This will allow us to test the discretization of the hyperbolic term and the treatment of the kinetic-coupling term. The numerical experiments were carried out assuming an exponential initial distribution, i.e., f*(V˜ ,0) ) f /0(V˜ ,0) ) e-V˜ , and using an uniform grid in r. In the simulations done with the 0-1-S model, κ was set equal to 1, which corresponds to a steady-state average number of radicals per particle, nj ) 1/3 (nj ) λ1). In the 0-1-2 simulations, γ was set equal to 2, whereas κ was calculated to obtain the same value of nj () λ1 + 2λ2). Finally, two values of R were considered, to evince the effect of stochastic broadening. Of course, such choices are arbitrary. Table 6 summarizes the values of parameters used in the various simulations as well as the errors of the solutions. Additionally, Figure 2 depicts the time-evolution of the overall distribution for two representative simulations. The plots relative to the individual distributions (not shown) are very similar to the ones presented here, since in this case f /m ≈ λmf*. With respect to the 0-1-S model, Figure 2 and Table 6 show that the agreement between numerical and analytical solutions is very good, even for sharp particle size distributions. We also see that the solutions are completely free from spurious oscillations (dispersion), contrary to what would have happened if high-order linear FD schemes had been employed. Note that

Figure 2. (case 1) 0-1-S model. Comparison between numerical (O) and analytical (s) solutions: (a) R ) 1, M ) 100; (b) R ) 10, M ) 200.

the breadth of the distributions narrows with time, a consequence of the fact that the radial growth rate, 1/(3r˜2), decreases with particle size. This effect is still more significant at higher values of R, because the extent of stochastic broadening diminishes. Judging from the errors in the moments, the quality of the solutions is equally good in the case of the 0-1-2 model. Note that the deviations are similar to those observed with the 0-1-S model at equal values of R and M. Besides, as γ f ∞, the solutions obtained with the 0-1-2 model converge to the solutions of the 0-1-S model, thereby demonstrating the consistency of the algorithm.

Ind. Eng. Chem. Res., Vol. 46, No. 2, 2007 649 Table 7. Results for Case 2 τ 10

50

100

0-1-S: R ) 1, κ ) 1, Ω ) 1, M ) 100 0.04 0.04 L2 error f /m (%) error λm (%) 0.00 0.00 0.04 0.01 error µm (%)

0.05 0.00 0.01

0-1-S: R ) 10, κ ) 1, Ω ) 1, M ) 100 0.34 0.95 L2 error f /m (%) error λm (%) 0.00 0.00 error µm (%) 0.15 0.14

1.36 0.00 0.12

0-1-2: R ) 10, κ ) 2.11, γ ) 2, Ω ) 1, M ) 100 0.30 0.88 1.28 L2 error f /m (%) error λm (%) 0.00 0.00 0.00 error µm (%) 0.14 0.14 0.12

4.2. Case 2. Although this case is physically more complex than Case 1, the initial conditions are simpler, which makes it possible to obtain f /m solutions for the 0-1-2 PBEs as well. This case will thus permit a more complete evaluation of the quality of the solutions obtained for the 0-1-2 model, as well as a check of the treatment of the nucleation term. The values of the parameters employed in the numerical simulations and the corresponding errors are listed in Table 7. Additionally, two examples of the quality of the solutions are presented in Figure 3. These results demonstrate that the numerical solutions are in excellent agreement with the analytical solutions for both models. In particular, the steep moving front is accurately reproduced: there is no dispersion and little dissipation is observed, even for very sharp fronts (Figure 3b). Note also, as pointed out in Section 1, the ability of the method to treat singular source terms in a straightforward manner, without the need to approximate the Dirac delta function by some other function. This ability translates into an exact preservation of the total number of particles (λ), which is in turn reflected in precise values of the individual moments λm. 4.3. Case 3. This case allows us to test the technique proposed to discretize the aggregation terms. The simulations presented here were done taking for initial condition a seed composed of particles with one radical and with an exponential size distribution, i.e., f*(V˜ ,0) ) f /1(V˜ ,0) ) e-V˜ . Also, a geometrical mesh in r was used, since this type of grid is the most appropriate for aggregation-dominated problems. A comparison of numerical and analytical results is depicted in Figures 4 and 5. For readability, the plots represent normalized distributions, f /m/λm. A summary of the errors can be found in Table 8. As demonstrated by Figures 4 and 5 and by the low L2-errors, the agreement between the two solutions is very good throughout the entire domain. Additionally, we can see that the moments λm and µm are accurately predicted, thus, corroborating the internal consistency of the generalized fixed pivot technique described in Section 3.2 with respect to these two properties. 4.4. Case 4. The purpose of this last idealized situation is to demonstrate the capacity of the method to deal with the simultaneous occurrence of the phenomena analyzed independently in Cases 1-3. Given the impossibility of deriving analytical solutions for the distributions, the comparisons must be made exclusively in terms of moment values. Figure 6 shows the simulation results for the 0-1-2 model. The curves obtained with the 0-1-S model are very similar, since the parameters κ and γ were chosen to give the same steady-state nj value. Initially, the number of particles increases (as indicated by the area below the curves), but progressively, an equilibrium between nucleation and aggregation is established

Figure 3. (case 2) Comparison between numerical (O) and analytical (s) solutions: (a) 0-1-S model, R ) 1; (b) 0-1-2 model, R ) 10.

and the particle number stabilizes. As can be seen in Table 9, the deviations in λm are insignificant and the errors in µm are comparable to the ones registered in Case 2. Judging from the good predictions of the moments, we may infer (indirectly but in all probability) that the distributions computed numerically are themselves correct and in good accord with the true solutions. 4.5. Real Case Study. To demonstrate the value of the technique described in this paper, we present some results of simulations relative to the emulsion polymerization of vinyl chloride. These simulations were carried out with a first

650


Figure 4. (case 3) 0-1-S model. Comparison between numerical (O) and analytical (s) solutions.

principles population balance model based on the 0-1-2 PBEs. This is actually the first time that the 0-1-2 model is used to compute the evolution of PSD; its interest and advantages will be discussed thoroughly in a subsequent publication. The example given here refers to an ab initio polymerization conducted in a batch reactor. Particle nucleation is assumed to occur by micellar and homogeneous nucleation, and particle coagulation is modeled following the work of Arau´jo et al.26 The evolution of the PSD and the number of particles is depicted in Figure 7. As a major feature, we observe that the particle number initially rises but then diminishes because the small

Figure 5. (case 3) 0-1-2 model. Comparison between numerical (O) and analytical (s) solutions. Results for f /1 are identical to those shown in Figure 4b.

precursor particles so formed are insufficiently stable and coagulate. An equilibrium number of particles is reached when the particles acquire enough stability. These simulations were done with 200 cells and took about 3.5 min on a 2.4-GHz PC. Approximately 90% of the CPU time is consumed in the aggregation terms. Of course, if automatic time-step control was implemented, computational requirements would be much inferior. In any case, CPU times of this magnitude are perfectly affordable even with a desktop PC,


Figure 6. (case 4) 0-1-2 model, numerical solutions. Results for the 0-1-S model are very similar. Table 8. Results for Case 3 τ 10

50

100


0-1-S: y10 ) 1, M ) 100 0.20 0.25 0.00 0.00 0.00 0.00

0.26 0.00 0.00


0-1-2: y10 ) 1, M ) 100 0.25 0.25 0.00 0.00 0.00 0.00

0.26 0.00 0.00

Table 9. Results for Case 4 τ 2

5

0-1-S: R ) 1, κ ) 1, Ω ) 1, Λ ) 0.2, M ) 100 error λm (%) 0.00 0.00 error µm (%) 0.30 0.09

10 0.00 0.05

0-1-2: R ) 1, κ ) 2.11, γ ) 2, Ω ) 1, Λ ) 0.2, M ) 100 error λm (%) 0.00 0.00 0.00 error µm (%) 0.20 0.07 0.04

therefore demonstrating that the proposed numerical method is a viable option for simulating emulsion polymerization processes. 5. Conclusions In this work, we have developed a method to numerically solve the coupled systems of integro-hyperbolic partial differential equations that appear in the 0-1 and 0-1-2 population balance models of emulsion polymerization. The advective term was discretized by means of a high-resolution FV scheme, which avoids unphysical oscillations and minimizes numerical diffusion. The fixed pivot technique of Kumar and Ramkrishna was extended to systems of PBEs and used to approximate the aggregation integrals in a way that ensures internal consistency with regard to particle number and mass. Implicit-explicit methods were employed to deal with the different natures of

Figure 7. Simulation of the synthesis of a PVC latex: (a) evolution of the PSD; (b) evolution of the number of latex particles.

the source terms. Additionally, new analytical solutions (distributions and/or their moments) were determined. Four cases were considered to evaluate the performance of the algorithm: (i) seed growth, (ii) nucleation and growth, (iii) aggregation, and (iv) simultaneous nucleation, growth, and aggregation. In all cases, numerical and analytical solutions were found to agree rather well. Additionally, simulations for a real case study demonstrated the ability of the technique to deal with physically realistic parameters and functions and, thus, its value for practical applications. A natural improvement to the method would be the inclusion of an adaptive time-step strategy. As mentioned in Section 3.3, this only depends on having good error estimates, which should be available soon (this is a very active field of research).

652


Finally, it is worth stressing that although the method was developed having in mind emulsion polymerization systems, it is actually quite general and can therefore be applied to any other processes governed by PBEs with similar advective and source terms.

Table B1. PBE Solutions for the 0-1-S Modela case 1b

Acknowledgment The authors express their gratitude to W. Hundsdorfer (CWI) for useful discussions on IMEX methods and to A. Mendes and F. Magalha˜es (University of Porto) for the careful review of the manuscript. H.M.V. thanks the Portuguese Science and Technology Foundation for financial support (grant SFRH/BD/ 10513/2002). Appendix A: Derivation of Eqs 13, 14, and 16 A.1. Equations 13 and 14. Dividing both sides of eq 12 by ∆ri and making use of eq 11a, we obtain

1

(1 - 1/2δmq)

∆ri

∫R dV ηi(V)∫VV-V i

nuc

V′,V′)

dV′ β(V -

nuc

M

F /0(s,τ) )

F /0(s,0) τ exp -(R(2 + κ) + s - b) 2b 2 [(b + Rκ + s) + (b - Rκ - s)e-bτ]

F /1(s,τ) )

RF /0(s,0) τ exp -(R(2 + κ) + s - b) (1 - e-bτ) b 2

3

(A1)

)

1

∑ ∫R dV ηi(V)∫V

(1 - 1/2δmq) ∆ri j,k)1

V-Vnuc

i

F /1(V˜ ,s) )

Ω exp s

F*(s,τ) )

4F*(s,0) (τ + 2)[τ + 2 - τF*(s,0)]

- V′ -

1

4F /1(s,0) [τ + 2 - τF*(s,0)][τ + 2 + 2τF /1(s,0) - τF*(s,0)]

4 2V˜ exp τ+2 (τ + 2)2

(

)

(

y1,0τ f /1(V˜ ,τ) ) f *(V˜ ,τ) 1/2 + (y1,0 - 1/2) exp -2V˜ τ+2

- Vk) (A2)

Applying the shifting property of the Dirac delta function, we find that the only nonzero terms of the sum are those for which Vˆ ) (Vj + Vk) ∈ Ri and, thus,

B h ri )

)]

R2(1 + κ) -R(1 + κ) - s V˜ s+R

[

dV′ β(V -

Vj)N qk δ(V′

[(

IC: f *(V˜ ,0) ) e-V˜ , f /1(V˜ ,0) ) y1,0f*(V˜ ,0)

nuc

V′,V′)Nmj δ(V

]

R(1 + κ) / F (V˜ ,s) s+R 1

f*(V˜ ,τ) )

B h ri

[

F /0(V˜ ,s) )

F /1(s,τ) )

We now interchange the order of integration and summation to get M

]

BC: f /1(0,τ) ) Ω IC: f /m(V˜ ,0) ) 0

M

Nmj δ(V - V′ - Vj)∑ Nqk δ(V′ - Vk) ∑ j)1 k)1

[

b ) x[R(2 + κ) + s]2 - 4Rs 2c

B h ri )

BC: f /m(0,τ) ) 0 / IC: f m*0 (V˜ ,0) ) 0

)]

a F / ) L [f / ]; Laplace transforms may either be with respect to V ˜ m m (cases 1 and 3) or with respect to τ (case 2). b Equivalent solutions, but in terms of the Fourier transform, have been obtained by Lichti et al.17 c Equivalent solutions have been obtained by Lichti et al.17

Table B2. PBE Solutions for the 0-1-2 Model

∑

(1 - 1/2δmq) ηi(Vˆ )β(Vj,Vk)Nmj Nqk ∆ri j,k∈Qi

Qi ) {{j,k}: 1 e k, j e i; Vi-1 < Vˆ < Vi+1}

(A3)

For m * q, this equation results in eq 13. On the other hand, for m ) q, we obtain a symmetrical sum that can be further simplified by removing the repeated combinations, leading to

B h ri )

1

∑

∆ri j,k∈Qi

(1 - 1/2δjk)ηi(Vˆ )β(Vj,Vk)Nmj Nmk

Qi ) {{j,k}: 1 e k e i; k e j e i; Vi-1 < Vˆ < Vi+1}

(A4)

A.2. Equation 16. A similar reasoning is used for the death term. We divide both sides of eq 15 by ∆ri and make use of eq 11 to obtain

D h ri ) -

1

∞ dV Nmi δ(V - Vi)∫V ∫ C ∆r i

M

dV′ β(V,V′)

nuc

∑ Nkδ(V′ -

k)1

i

Vk) (A5) and then interchange the order of integration and summation to get

D h ri ) -

a

Solutions could only be derived for the initial condition y1,0 ) 1.

Finally, by using the shifting property of the Dirac delta function, we arrive at

M

1

∞ dV N mi δ(V - Vi)∑ ∫V ∫ C ∆r i

i

k)1

nuc

dV′ β(V,V′)Nkδ(V′ Vk) (A6)

D h ri ) -

1 ∆ri

M

Nmi

∑β(Vi,Vk)Nk

k)1

(A7)


B4. Here, the symbol “≈” denotes steady-state or long-time limit solutions (valid for the τ values used in this work).

Table B3. Moment Solutions for the 0-1-S Model case λ ) 1

1

λ0 )

Nomenclature

1 + κ + exp[-R(2 + κ)τ] 2+κ

µ)1+

[

1 1 τ+ (exp[-R(2 + κ)τ] - 1) 2+κ R(2 + κ)

[

1+κ 2 τ + (2 + κ) R(2 + κ) (2 + κ)2 λ ) Ωτ µ0 ≈

2

λ0 ) µ≈

[

λ)

[R(2 + κ)τ + e-R(2+κ)τ - 1]

(1 + κ) Ω τ2 (1 + κ)τ + (2 + κ) 2 R(2 + κ) R2(2 + κ)2

µ0 ≈ 3

Ω(1 + κ) R(2 + κ)2

]

]

[

]

Ω(1 + κ) τ2 (1 + 2κ) κτ + (2 + κ)2 2 R(2 + κ) R2(2 + κ)2

]

2 2+τ

y1,0 1 + y1,0τ µ ) 1

λ1 )

µ1 )

λ)

4

y1,0(2 + 2τ + y1,0τ2) 2(1 + y1,0τ)2

x2ΩΛ tanh(xΩΛ2 τ)

Table B4. Moment Solutions for the 0-1-2 Model case 1

2

λ ) 1 λ1 ≈

1 + 2(κ + γ) (3 + 2κ + 2γ)κ + 4γ + 2

λ2 ≈

1 (3 + 2κ + 2γ)κ + 4γ + 2

λ ) Ωτ λ1 ≈

[

Ω (1 + 2κ + 2γ)τ + b

]

4κ3 + 4(1 + 2γ)κ2 + (1 + 4γ2 + 8γ)κ + 4γ(γ + 1) + 1 bR λ2 ≈

[

]

2κ(γ + κ) + 2γ -1 Ω τ+ b bR

where: b ) (3 + 2κ + 2γ)κ + 4γ + 2 3

µ2 ) 4

Greek Letters

λ, λ1, µ, and µ1: see Table B3 (y1,0 ) 1) λ2 )

x2 -1 x2 (1 + τ) x2 2(τ + 1) (1 + τ) +1

4τ(2 + τ)(1 + τ)

x2

x2

+ x2(2 + τ(2 + τ))((1 + τ)2

4[1 + τ + (1 + τ)

a(f,x,y) ) vector of aggregation terms (part m-3 [x]-1 s-1) B ) generic birth term (part m-3 [x]-1 s-1) c ) pseudo-first-order rate coefficient for termination (s-1) D ) generic death term (part m-3 [x]-1 s-1) e ) vector of terms to be integrated explicitly (part m-3 m-1 s-1) f(x,t) ) vector number density function (part m-3 [x]-1) f (x,t) ) overall number density function (part m-3 [x]-1) F ) Laplace transform of f fm(x,t) ) number density function for particles having radicals of type m (part m-3 [x]-1) h ) flux function (part m-3 s-1) kdes ) desorption frequency (s-1) kdM ) rate coefficient for desorption of monomeric radicals from particles (s-1) kfM ) rate coefficient for transfer to monomer (m3 mol-1 s-1) kp′ ) rate coefficient for propagation of a monomeric radical (m3 mol-1 s-1) K ) rate coefficient of volume growth (m3 s-1) M ) number of finite volumes [M]p ) concentration of monomer in a particle (mol m-3) nj ) average number of radicals per particle n ) nucleation vector Ni ) total number of particles in cell i (part m-3) Nmi ) number of particles of type m in cell i (part m-3) r ) unswollen radius of a particle (m) ri ) center of the ith cell (m) ∆ri ) size of the ith cell (m) r3 ) vector of radial growth rates (m s-1) R ) kinetic coupling matrix (s-1) Rnuc ) rate of particle nucleation (part m-3 s-1) s(f,x,y)) source function (part m-3 [x]-1 s-1) t ) time (s) V ) unswollen volume of a particle (m3) Vˆ ) unswollen volume of a new particle formed by aggregation (m3) Vi ) unswollen volume corresponding to ri (m3) ∆Vi ) size of the ith cell (m3) v3 ) vector of volume growth rates (m3 s-1) Vw ) volume of the aqueous phase (m3) y ) continuous-phase vector ym,0 ) initial number fraction of particles with m radicals

- 1)

1+x2 2

]

see Table B3

Appendix B: Analytical Solutions Solutions to the PBEs of Table 3, obtained by the method of Laplace transforms,6 are summarized in Tables B1 and B2. For the cases where the moment equations of Table 4 could be integrated analytically, solutions are given in Tables B3 and

R ) ratio of entry to growth β ) coagulation rate coefficient (m3 part-1 s-1) γ ) ratio of termination to entry δ(x) ) Dirac delta function ([x]-1) δij ) Kronecker delta ηi ) particle fraction assigned to cell i κ ) ratio of desorption to entry λ ) total number of particles λm ) number of particles with m radicals Λ ) ratio of coagulation to growth µ ) total volume of particles µm ) volume of particles with m radicals θ ) residence time (s-1) F ) total entry frequency of radicals into particles (s-1) FE ) entry frequency of monomeric radicals into particles (s-1)

654


FI ) entry frequency of initiator-derived radicals into particles (s-1) τ ) dimensionless time Ω ) ratio of nucleation to growth Subscripts, Superscripts, and Accents in ) reactor inlet max ) upper bound M ) monomeric radical nuc ) nucleated particles out ) reactor outlet P ) polymeric radical r ) relative to radius v ) relative to volume - ) cell average value ∼ ) dimensionless coordinate * ) dimensionless number density function Literature Cited (1) Vale, H. M.; McKenna, T. F. Modeling particle size distribution in emulsion polymerization reactors. Prog. Polym. Sci. 2005, 30, 1019-1048. (2) Gilbert, R. G. Emulsion polymerization. A mechanistic approach; Academic Press: San Diego, 1995. (3) Lichti, G.; Gilbert, R. G.; Napper, D. H. Molecular weight distribution in emulsion polymerizations. J. Polym. Sci., Part A: Polym. Chem. 1980, 18 (4), 1297-1323. (4) Clay, P. A.; Gilbert, R. G. Molecular weight distributions in freeradical polymerizations. 1. Model development and implications for data interpretation. Macromolecules 1995, 28, 552-569. (5) Prescott, S. W.; Ballard, M. J.; Gilbert, R. G. Average termination rate coefficients in emulsion polymerization: effect of compartmentalization on free-radical lifetimes. J. Polym. Sci., Part A: Polym. Chem. 2005, 43, 1076-1089. (6) Ramkrishna, D. Population balances. Theory and applications to particulate systems in engineering; Academic Press: San Diego, 2000. (7) Immanuel, C. D.; Doyle, F. J., III Computationally efficient solution of population balance models incorporating nucleation, growth and coagulation: application to emulsion polymerization. Chem. Eng. Sci. 2003, 58, 3681-3698. (8) Attarakih, M. M.; Bart, H.-J.; Faqir, N. M. Numerical solution of the spatially distributed population balance equation describing the hydrodynamics of interacting liquid-liquid dispersions. Chem. Eng. Sci. 2004, 59, 2567-2592. (9) Alexopoulos, A. H.; Roussos, A. I.; Kiparissides, C. Part I: Dynamic evolution of the particle size distribution in particulate processes undergoing combined particle growth and aggregation. Chem. Eng. Sci. 2004, 59, 57515769. (10) Alexopoulos, A. H.; Kiparissides, C. Part II: Dynamic evolution of the particle size distribution in particulate processes undergoing simultaneous particle nucleation, growth and aggregation. Chem. Eng. Sci. 2005, 60, 4157-4169.

(11) LeVeque, R. J. Finite Volume methods for hyperbolic problems; Cambridge University Press: Cambridge, 2002. (12) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization-I. A fixed pivot technique. Chem. Eng. Sci. 1996, 51 (8), 1311-1332. (13) Lim, Y. I.; Le Lann, J. M.; Meyer, X. M.; Joulia, X.; Lee, G.; Yoon, E. S. On the solution of population balance equations (PBE) with accurate front tracking methods in practical crystallization processes. Chem. Eng. Sci. 2002, 57, 3715-3732. (14) Roussos, A. I.; Alexopoulos, A. H.; Kiparissides, C. Part III: Dynamic evolution of the particle size distribution in batch and continuous particulate processes: A Galerkin on finite elements approach. Chem. Eng. Sci. 2005, 60, 6998-7010. (15) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization-III. Nucleation, growth and aggregation of particles. Chem. Eng. Sci. 1997, 52 (24), 4659-4679. (16) Coen, E. M.; Gilbert, R. G.; Morrison, B. R.; Leube, H.; Peach, S. Modelling particle size distributions and secondary particle formation in emulsion polymerization. Polymer 1998, 39 (26), 7099-7112. (17) Lichti, G.; Gilbert, R. G.; Napper, D. H. The growth of polymer colloids. J. Polym. Sci., Polym. Chem. Ed. 1977, 15, 1957-1971. (18) Hundsdorfer, W.; Verwer, J. G. Numerical solution of timedependent adVection-diffusion-reaction equations; Springer-Verlag: Berlin, 2003. (19) Liu, X.-D.; Osher, S.; Chan, T. Weighted essentially nonoscillatory schemes. J. Comput. Phys. 1994, 115, 200-212. (20) Shu, C.-W. Essentially nonoscillatory and weighted essentially nonoscillatory schemes for hyperbolic conserVation laws; NASA/CR-97206253, ICASE Report No. 97-65. 1997. (21) Ascher, U. M.; Ruuth, S. J.; Wetton, B. T. R. Implicit-explicit methods for time-dependent PDEs. SIAM J. Numer. Anal. 1995, 32, 797823. (22) Wang, D. Variable step-size implicit-explicit linear multistep methods for time-dependent PDEs. MSc Thesis, Simon Fraser University, Burnaby, 2005. (23) Feeney, P. J.; Napper, D. H.; Gilbert, R. G. The determinants of latex monodispersity in emulsion polymerization. J. Colloid Interface Sci. 1987, 118, 493-505. (24) Feeney, P. J.; Napper, D. H.; Gilbert, R. G. Coagulative nucleation and particle size distributions in emulsion polymerization. Macromolecules 1984, 17 (12), 2520-2529. (25) Valko´, P. P.; Abate, J. Comparison of sequence accelerators for the Gaver method of numerical Laplace transform inversion. Comput. Math. Appl. 2004, 48, 629-636. (26) Arau´jo, P. H.; de la Cal, J. C.; Asu´a, J. M.; Pinto, J. C. Modeling particle size distribution (PSD) in emulsion copolymerization reactions in a continuous loop reactor. Macromol. Theory Simul. 2001, 10, 769-779.

ReceiVed for reView July 17, 2006 ReVised manuscript receiVed October 4, 2006 Accepted October 4, 2006 IE060928L

Solution of Population Balance Equations for Emulsion Polymerization

Recommend Documents