Prediction of Multimodal Distributions in Breakage Processes

Ind. Eng. Chem. Res. 2005, 44, 2649-2658

2649

Prediction of Multimodal Distributions in Breakage Processes P. Canu* Dipartimento di Principi e Impianti di Ingegneria Chimica, Universita` di Padova, Via Marzolo 9, 35131 Padova, Italy

The prediction of multimodal distributions through the population balance equation (PBE) in breakage processes is addressed. A numerical approach based on physical arguments is suggested and applied. The time-dependent unknown distribution is approximated by a weighted sum of β-subdistributions. The parameters of the β-PDFs and their weights are treated as unknown functions of time. The method allows for the tracking of the evolution of multimodal distributions through few unknown functions, determined by ordinary differential equations. The method is applied here to PBEs for pure breakage processes, although it can be further developed to include growth and coalescence. Introduction The population balance equation (PBE) is a paradigm in the quantitative description of many processes. It is a balance equation applied to a quantity that depends on several variables, through a distribution. It fits perfectly (no alternatives) the case of populations of many entities that are treated together, still having individual, clearly identifiable features, such as particles, bubbles, droplets, macromolecules, etc. Polymer reaction engineering is based on this equation, as the polymer characteristics (chain length, composition, branching degree) are intrinsically different among molecules that are all part of the same polymer sample.1 The PBE specifically accounts for a distinct feature of the populations, i.e., internal rearrangement processes. These are transformations that affect the so-called2 internal variables. Size is a typical internal coordinate, and this work addresses the processes of size reduction, broadly called “breakage” processes. Many applications can be classified as breakage processes, including polymer degradation3 and particle grinding.4 Refrences to other applications can be found in Ramkrishna’s book2 and many works5,6 discussing the solution of the PBE in a general context. The internal rearrangements inside a population can result in a variety of behaviors. Quite often, broadening and displacement of the distribution during its evolution can be observed. A distribution broadens when its members transform in many different way, depending on their size. When either the process evolves significantly from its initial state or the internal rearrangements move or create population members very far from their original positions along the internal coordinate, then multimodality can arise. Multimodality is the identification of different subpopulations with more or less distinct features. Typically, multiple local maxima emerge in the distribution. Collective properties of a multimodal population can be dramatically different. The additional details in multimodal distributions introduce some difficulties in the PBE solution methods. The simplest approach, based on the discretization of the distribution, requires far more size classes to reproduce multimodal behavior. On the other hand, * Tel.: +39 049 8275463. Fax: +39 049 8275461. E-mail: [email protected].

methods of solution based on continuous approximations must be able to describe functions with multiple maxima, which is not the case of standard probability distributions, such as the log-normal, Γ, and Weibull distributions, among others. Here, we introduce a new method within this approach, based on the observation that a multimodal distribution can be seen as a combination of subdistributions. Typically, analytical probability density functions (PDFs) are used as continuous approximations. Actual distributions can be well fitted, provided some sort of “perturbation” of the analytical PDF in terms of problem-specific polynomials, usually developed from the same analytical PDF, by orthogonalization.5,26 In our experience,7 such a method can be very efficient as long as the simulated distribution is monomodal and does not change too much with respect to the base PDF, which is not the case when multimodality develops. In this work, we suggest a method to use the same approach, based on a continuous PDFbased approximation, that can follow the rise and evolution of multimodality. We suggest the selection of an analytical PDF flexible enough to fit dramatic modifications in the actual PDF, simply by adjusting its parameters, and the combination of several such PDFs to describe very broad or multimodal distributions. It was already observed8,9 that multimodal distributions can be approximated by a suitable combinations of standard (log-normal, normal, and also β) PDFs. To our knowledge, only a few researchers10,11 have attempted to use a weighted sum of PDFs in simulation. They either impose, a priori, the time variation of the parameters or keep some of them constant, based on physical arguments. Accordingly, the evolution of the distribution is simply mimicked, not actually predicted from a conservation principle, which is the breakage PBE. In the polymerization literature ,an approach comparable to the one proposed here has been suggested.25 It is based on what is called there a “numerical fractionation” of the overall distribution in a weighted number of classes. Three moments (zeroth, first, and second) are then calculated for each class, using standard closure equations26 to calculate the third moment; from the three moments of each class, the PDF is reconstructed using an analytical PDF obtained for simpler configurations (Flory’s25 or Wesslau’s27 distribu-

10.1021/ie049652c CCC: $30.25 © 2005 American Chemical Society Published on Web 12/23/2004

2650

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005

tion); the complete PDF is then a sum of the PDFs of each class, weighted by the first moment. The technique has been successfully used to model the degradation of polystyrene,28 which can involve bimodal distributions. Differences from the method presented here will be detailed later, once our technique has been introduced. Here, we develop and test a method that uses a combination of standard, flexible probability distribution functions (PDFs) as tentative solution of the breakage PBE. The actual solution is achieved by reformulating the PBE as a set of ordinary differential equations that predict the time variation of the PDF parameters. A reduction of the number of equations required to predict the dynamics of a PBE is particularly appealing in applications where the external variables can vary significantly. This usually calls for a coupling of PBE and computational fluid dynamics (CFD) solvers, where a detailed description at two different scales (fluid flow and local PBE) definitely requires dedicated methods. PBE for Breakage In this section, we define the single form of the PBE, the solution of which will be illustrated in the following sections, under different circumstances and with different methods. It is worth anticipating that the equation and its solution methods have been written and developed for application to size-reduction processes of solid particles (milling, crushing, etc.). In this case, the internal coordinate (i.e., the size) corresponds to particle volume or mass. However, the same equation and solution methods can be equally well applied to other physical systems (bubbles, drops, polymers), even when the “size” is measured differently, such as in terms of chain length or composition in macromolecules. The PBE for a spatially homogeneous breakage process, without convective fluxes (closed system), can be synthetically written as2

∂f1(V,t) ) h(V,t) ) B(V,t) - D(V,t) ∂t

(1)

where f1(V,t) is the instantaneous number density distribution and B and D are “birth” and “death” functions, respectively. Accordingly, f1(V,t) dV measures the number of particles of size between V and V + dV at time t. The birth and death functions, B and D, account for the rate of acquisition and loss of particles in or from size V. They can be written as

B(V,t) )

∫V∞D(V′,t) ν(V) P(V|V′) dV′

D(V,t) ) b(V,t) f1(V,t)

(2)

Note that the rate of increase of the number of particles of size V, B(V,t), must follow the loss of particles from larger sizes, D(V′>V,t), multiplied by the probability P(V|V′) of yielding an average number ν(V) of particles (daughters) in the size V after breakage of a particle of size V′. The birth function B(V,t) involves an integral so that the PBE (eq 1) becomes an integropartial-differential equation in the unknown function f1(V,t). Solution of the PBE provides the evolution of the distribution f1(V,t) at any time. The initial distribution, f1(V,0); the daughter distribution function, ν(V)P(V|V′); and the breakage rate, b(V,t); collectively referred to as

breakage functions,2 must be specified to proceed with the solution, together with a suitable numerical method. The identification of the two functions ν(V) P(V|V′) and b(V,t) is largely dependent on the physics of the process. They contain the actual mechanism of breakage at the particle scale. The effort to develop efficient techniques for the solution of the PBE is mostly motivated by the inverse problem, i.e., the identification of the physics of the breakage [then expressions for ν(V) P(V|V′) and b(V,t)] from the fundamental modeling of actual, macroscopic data of f1(V,t). In the following discussion, rather general expression for ν(V) P(V|V′) and b(V,t) will be used. The breakage rate is usually assumed to increase with the size raised to some power R

b(V,t) ) kVR

(3)

where the parameter k and R can vary with the process conditions. They will be considered constant in the following discussion so that the breakage rate turns out to be time-invariant, b(V). It is important to highlight that the function b(V,t) appears in both functions B and D. Consequently, the constant k multiplies all of the rhs of the PBE so that it turns out to be a “kinetic” constant, that determines the overall rate of the breakage process. The form of ν(V) P(V|V′) and the size dependency of b(V)govern the rearrangements of particles inside the distribution, but the constant k specifies the extent of the evolution of the distribution from the initial configuration f1(V,0), at a given time. The daughter distribution function, ν(V) P(V|V′), is perhaps the most important part of the breakage equation, where the fragmentation mechanism is actually translated into the balance equation. A discussion of several breakage equations is available in the literature.17 A useful general form that can explicitly or approximately represent nearly all known breakage functions has been derived by Diemer and Olson.22 So far, we assume that ν(V) P(V|V′) has a self-similar behavior2 so that it can be written as

ν(V) P(V|V′) )

θ(z) V′

z≡

V V′

(4)

Valid daughter distributions must satisfy the normalization conditions and the conservation of mass after breakage.2 The breakage rate of eq 3 and daughter distribution functions of eq 4 above do not depend on f(x,t), which is a very useful feature for some methods of solution of the PBE. According to the general expressions of the breakage functions given above, the breakage PBE (eq 1) can be written as

∂f1(V,t) ) k[ ∂t

∫V∞(V′)R-1 f1(V′,t) θ(V/V′) dV′ - VR f1(V,t)]

(5)

A nice feature of pure breakage processes is the possibility of uniquely identifying a maximum size VM, which is the largest particle initially present. During the process, only particles smaller than or equal to VM can be observed. Accordingly, the particle size can be scaled to a 0-1 range defining the dimensionless variable x, and the PBE can be reformulated as

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2651

∂f1(xVM,t) ) k′[ ∂t

∫x1(x′)R-1f1(x′VM,t) θ(x/x′) dx′ xR f1(xV,t)]

x≡

V (6) VM

where k′ ) kVMR. In addition, it is worth transforming the number density f1 to a normalized distribution function, f, expressed in terms of mass density. The first requirement is useful for the method that will be developed later, which makes use of standard probability density functions (PDFs). The second requirement is partly motivated by the first (the total mass is the only conserved quantity), partly because classification of particles by mass is the customary formulation for experimental particle size distribution data in the milling literature, because of the measurement techniques used there. The connection between f1 and f is simplyf(x,t) ≡ Cxf1(xVM,t), where C is a constant ()FVM/Mtot) that will simplify. Accordingly

∫xx f(x,t) dx ) 2

1

M1-2 Mtot

(7)

measures the mass fraction between size x1 and x2, and consequently

∫01f(x,t) dx ) 1

(8)

The PBE can then be written in terms of the mass density distribution f

∫x1(x′)R-2f(x′,t) θ(x/x′) dx′ - xR f(x,t)]

∂f(x,t) ) k′[x ∂t

(9)

which is the form used throughout in the following discussion. The PBEs for pure breakage and number density can be analytically solved for some simplified expression of the breakage functions23 or in the form of similarity solutions.6,18,24 Unfortunately, they are typically monomodal and sometimes simply asymptotic solutions. The method developed in the following section is dedicated to multimodal distributions, which are usually transient configurations of practical interest. PDFs-Based Solution of the Breakage Equation In this section, we introduce a novel method of solution of the PBE, based on physical considerations about the expected solution. The interest is specifically in multimodal solutions, not excluding monomodal cases. Both can be observed in the same process, at different times. The method is based on assumptions about the explicit analytical form of the expected solution. The approximating function adjusts to the evolution dictated by the PBE by means of parameters that vary in time. The approximating solution, g(x,t,λ), must be flexible enough, through its parameters λ, to describe very different observable size distributions. To describe multimodal distributions, we assume that the function g(x,t,λ) approximating the real solution f(x,t) has the form

g(x,t,λ) )

∑i wi(t) gi(x,t,ηi)

(10)

where each component, gi(x,t,ηi), is a suitable subdistribution. It is likely, indeed, that a multimodal distri-

Figure 1. β-PDFs for different values of parameters p and q. Each PDF has been scaled by its maximum value for readability.

bution is the result of the coexistence of more than one “familiy” in the whole population. Note that λ includes both the parameters of each subdistribution, ηi, and the weights, w. The weighted sum of PDFs in eq 10 is not simply used to reconstruct the complete PDF as in the numerical fractionation technique,25 but it is directly used in the solution of the PBE, without going through its solution by moments.25,27 Parameters λ are determined by the solution of the complete PBE for g(x,t,λ), instead of being calculated from moments. Taking advantage of the closed size domain, limited to [0, 1], analytical probability density functions (PDFs) defined on [0, 1] can be used. Among others, the β-PDF has a very simple expression and can take extremely different shapes, as shown in Figure 1. Note that Gaussian shapes, possibly very narrow, skewed toward small or large sizes, and centered at widely different sizes, and exponential shapes are well approximated. These features of the β-PDF have been observed and used in other applications dealing with distributions, such as PDF modeling of turbulent reactive flows.20 Accordingly, we suggest that β-PDFs be employed to express the subdistributions gi(x,t)

gi(x,t,ηi) ) βi(x,pi(t),qi(t))

(11)

where the explicit form of βi is

βi(x,pi,qi) )

1 xpi-1(1 - x)qi-1 B(pi,qi)

B(p,q) )

Γ(p+q) Γ(p) Γ(q)

(12)

Note that, in a single β-PDF, the size dependency reduces to a simple power law, while the time dependency remains in the parameters pi(t) and qi(t). The approximation in eq 10 is not guaranteed to represent any function by simply increasing the number of terms, because β-PDFs are an incomplete set. However, its capability to satisfactorily describe experimental multimodal distributions will be demonstrated in the following section. An approximation of f(x,t) that is analytically explicit in the size, x, and specifically a β-PDF, allows for the

2652


closed-form evaluation of the integral that appears on the rhs of the breakage PBE. The power-law size dependency of the β-PDF matches very neatly the expression of b(x) ) k′xR and many forms of θ(z). Accordingly, extensive analytical preprocessing of the equations can be done, greatly simplifying the PBE and its solution. Let us illustrate this point in the simples case of just one β-PDF

g(x,t,η) ) β(x,t,η) )

1 xp-1(1 - x)q-1 B(p,q) η ) [p(t) q(t)] (13)

Replacing f(x,t) with β(x,p,q) in the breakage PBE (eq 9), we obtain

∂β(x,p,q) ) k′[x ∂t

∫x1(x′)R-2β(x′,p,q) θ(x/x′) dx′ R

x β(x,p,q)] (14) The time derivative of the β-PDF can be evaluated as

∂β(x,p,q) ∂β(x,p,q) dp ∂β(x,p,q) dq ) + ∂t ∂p dt ∂q dt

(15)

showing that the time variation of f(x,t) remains in the two unknown functions p(t) and q(t). The partial derivatives of the β-PDF with respect to the parameters have the following expressions

∂β(x,p,q) ) [log(x) - Ψ(p) + Ψ(p+q)]β(x,p, q) ∂p ∂β(x,p,q) ) [log (1 - x) - Ψ(q) + Ψ(p+q)]β(x,p, q) ∂q (16)

1 dΓ(z) Γ(z) dz

0 < γ e1

∫x1(x′)R-γ-1β(x′,p,q) dx′ -

xRβ(x,p,q)]

∫x1β(x′,p1,q) dx′ B(p2,q) β(x,p2,q)] (19)

where

βINC(x,p,q) ≡

∫0xβ(y,p,q) dy

(21)

where βINC(x) is the β-cumulated probability, also called the incomplete β-function. Eventually, by means of eqs 16, 17, and 19, we translate the single PBE (partial integro-differential equation) into a single ordinary differential equation (ODE) in two unknown functions, p(t) and q(t)

∂β(x,t) dp ∂β(x,t) dq k′ [(γ + 1)xγB(p1,q)[1 + ) ∂p dt ∂q dt B(p,q) βINC(x,p1,q)] - B(p2,q) β(x,p2,q)] (22) which then requires an additional equation to be uniquely solved. It must be observed that the manipulations in the rhs of eq 14 that led to the analytical evaluation of the integral are sensitive to the expressions for b(x) and θ(x/x′). Whereas there are no significant alternatives to the power-law form of b(x) (and indeed we already included it in the general PBE of eq 9), several daughter distribution functions have been suggested, on the basis of both physical arguments and empirical data. Analytical manipulations are possible for many of them. In the case of the simpler form of the binary β-function29

θ(z) ) 60z2(1 - z)2 ) 60(z2 - 2z3 + z4)

) 60x3 )

60

(23)

∫x1(x′)R-2β(x′,p,q) θ(x/x′) dx′

∫x1(x′)R-4β(x′,p,q)[1 - 2/x′ + 1/(x′)2] dx′ 2

∑ Qi B(p,q)i)0

(24)

where each term is similar to the one identified above for the power law θ(z)

(18)

rhs ) k′[(γ + 1)xγ

k′ [(γ + 1)xγB(p1,q) B(p,q)

∫x1β(y,p,q) dy ) 1 - βINC(x,p1,q)

x

which is frequently used for particle grinding.19 Interestingly, it simply introduces and additional exponent to the size, so that the rhs becomes

)

(20)

and the integral of a β-PDF over the upper portion of its complete range is

(17)

is the digamma function, that can be easily evaluated numerically, the same as its parent gamma function, Γ. Depending on the form of the daughter distribution, θ(x/x′), the β-PDF might allow many simplifications also in the rhs of eq 14, thanks to the combinations of powers. We demonstrate the manipulations for the power-law daughter distribution22

θ(x/x′) ) θ(z) ) (γ + 1)zγ-1

p2 ) p + R

expansion of the product highlights a sum-of-powers structure, so that the integral on the rhs of eq 14 is split into three analytical terms

where

Ψ(z) )

p1 ) p + R - γ - 1

Qi ) (- 1)i

()

2 3+i x B(pi,q)[1 - βINC(x,pi,q)] i pi ) p + R - 4 - i (25)

Similarly, the generalized expression of the daughter distribution derived by Diemer and Olson22

θ(z) )

wi n i

∑i B(a ,c )za -1(1 - z)c -1 i

i

i

(26)

i

can be processed for positive integral values of the


parameters ci, for which

(-1)c-1-i( ∑ i i)0

c-1

za-1(1 - z)c-1 )

c-1

)

za+c-2-i (27)

is again a sum of powers. The more general case of a sum of NB β-PDFs in the approximation g(x,t,λ) is a linear combination of the analytical solutions above, thanks to the fact that f(x,t) [and hence g(x,t,λ)] appears linearly (or in linear operators) in the PBE. This would not be the case for breakage functions that include f(x,t), although such an alternative has never been suggested. We omit the transcription of the analogue of eq 22 for many β-PDFs and other daughter distributions, but we must observe that the weights wi(t) are additional unknown functions, that add to the parameters of each β-PDF, pi(t) and qi(t). In summary, we translated a single PBE into a single ODE and one constraint for the weights, Σiwi ) 1, i.e., two equations for 3NB unknown functions, pi(t), qi(t), and wi(t) for i ) 1, ..., NB. It is now clear that the method provides a dramatic analytical simplifications in the PBE, but leave us with at most two (if we differentiate the constraint on the weights wi) ODEs in many (3NB - 1) unknown functions. The latter problem can be solved by any method of weighted residuals. The method of moments nicely combines with the β-PDF, because it simply modifies the parameter p of any β-PDF as p + k, for any kth moment.12 Accordingly, equations for all of the moments can be formulated analytical for many breakage kernels. Moments must be translated back into the original parameters of the PDFs to provide a complete representation of the distribution. Eventually, a collocation method turns out to be simpler and allows for the direct manipulation of the whole distribution. It simply reduces to application of the ODE to a number, NC, of points xj in the range [0, 1]. In practice, a set of differential equations is generated from a single ODE written at several locations xi. When a power-law daughter distribution holds (eq 18) and a single β-PDF is being used in the approximant of eq 10, the single ODE obtained after preprocessing, eq 22, splits into NC ODEs

∂β(xi,t) dp ∂β(xi,t) dq k′ {(γ + + ) ∂p dt ∂q dt B(p,q) 1)xiγB(p1,q)[1 - βINC(xi,p1,q)] - B(p2,q) β(xi,p2,q)} i ) 1, ..., NC (28) when applied to NC collocation points xi. The NC ODEs can be written in compact matrix notation as

M(x,λ,t)

dλ ) N(x,λ,t) dt

(29)

where λ is the vector of NP unknown parameters [p(t) and q(t) for a single β-distribution; pi(t), qi(t), and wi(t) in the most general case]; matrix M contains the partial derivatives of β(x,t) [or g(x,t,λ) in the most general case] with respect to each parameter (in the columns), evaluated at each collocation point (along the row), resulting in a NC × NP size; and vector N contains the rhs of eq 22 evaluated at each collocation point, becoming an NC × 1 vector. Matrix M must be square and nonsingular to uniquely solve the linear system of eqs 29 and

Figure 2. Approximation of an experimental bimodal distribution with two β-PDFs. Data from ref 13. PVA ground in a shaker bead mill.

determine the time derivative of the parameters. This means NC ) NP ) 3NB - 1 distinct collocation points. Interestingly, NC can be higher than the number of parameters NP; system of eqs 29 is then overdetermined, and its solution is no longer unique, but a least-squares solution can be easily achieved. We will take advantage of this feature in the following calculations. We will show in the following section that multimodal distributions can be described very accurately with as few as two β-PDFs and five collocations points. Results Before we examine the results of solving the breakage PBE by the method described above, we need to confirm a basic assumption in it, i.e., the possibility that actual (multimodal) distributions can be accurately described by means of a finite weighted sum of β-PDFs. These PDFs are not expected to form a basis of functions, and therefore to guarantee, a priori, that any sort of function can be approximated by a sum of members of the basis, provided a sufficient number is taken. We need to provide evidence that a similar behavior is observed with β-PDFs as well, in actual circumstances. We tried to approximate experimental multimodal distributions by means of the approximant of eq 10. In practice, we have to determine the 3NB - 1 parameters λ in eq 10 that allow a fit to a set of experimental data much larger than the set of unknown parameters. This is done by a nonlinear least-squares technique applied to all of the experimental data available for each single size distribution. Figure 2 shows an experimental13 mass PDF of a polymer [poly(vinyl acetate), PVA] after milling. In this case, because of the nature of the material (aggregates), milling yields two clearly distinguished populations, fairly symmetric and approximately Gaussian in the logarithm of the size. It is a simple and clear example of a multimodal (actually, bimodal) distribution. Quite obviously, at least two β-PDFs are required to approximate these data. The approximation, as shown in Figure 2, turns out to be very good, also considering the experimental error. The first subpopulation, β1(x), describes the distribution of smaller particles and is almost zero elsewhere, and vice versa for the second, β2(x). Figure 3 shows a distribution in which multimodal

2654


Figure 3. Approximation of an experimental bimodal distribution with two β-PDFs. Data from ref 14. Al2O3‚3H2O ground in a stirred ball mill.

Figure 4. Approximation of an experimental bimodal distribution with three β-PDFs. Data from ref 14. Al2O3‚3H2O ground in a stirred ball mill.

features are less clearly shaped. Two β-PDFs provide a reasonably good approximation, although small discrepancies can be observed at the largest sizes (rightmost tail) and close to the maximum of the most abundant subpopulation, where the measured distribution cuts slightly lower than the approximation. It is interesting to observe how the two approximating function merge and skew to provide a weighted sum that closely matches the experimental data. Small inaccuracies can be reduced by increasing the number of subdistributions, as shown in Figure 4. Note that the physical basis of eq 10, i.e., suggesting an approximation function that contains subdistributions to model multimodal data, smoothly evolves into the use of β-PDFs as generic, flexible functions to provide a better fit. From Figure 4, one wonders whether the sample could actually be a mixture of three families of particles, each one with its own distribution, or whether that is simply a mathematical artifact. Figure 5 raises even more such a doubt: Does the middle distribution have a physical interpretation, or does the third β-PDF simply allow the mathematical solution to follow some experimental inaccuracy in the measurements?

Figure 5. Approximation of an experimental distribution with three β-PDFs. Data from ref 9. TiO2 wet ground in a stirred bead mill.

In any case, we believe that Figures 2-5 clearly demonstrate the flexibility of a weighted sum of β-PDFs and its impressive ability to describe actual multimodal distributions. Similar results have already been observed by others,8,9 also using combinations of lognormal PDFs. However, we now want to make a further step, from simply using many standard PDFs to fit experimental data, to using them in simulations through the pertinent PBE. Accordingly, we developed the method described above, and here, we show the results of its application. All of the calculation were checked against another numerical solution instead of analytical results, which are not available for calculations that involve transient multimodal distributions. The reference solution was obtained by discretization of f(x,t) (method of classes) and solution of the ODEs in each fraction, fi(xi,t). The integral in the birth term was solved by assuming constant fi(xi,t) in the size interval xi ( ∆x/2, without any correction for this approximation21 because we used a large number (N ) 200) of classes to guarantee accuracy. The discretized breakage PBE of eq 9 then reduces to

dfi dt

N

) k′[xi

xjR-2fj(t) θ(xi/xj)∆xj - xiRfi(t)] ∑ j)1

(23)

where the subscript i (or j) means that the corresponding function is evaluated at the ith (or jth) size. If the breakage function does not depend on f(x,t) or on time, as in eqs 3 and 4, the discretized PBE (eq 23) can be conveniently rearranged in vector notation as

df ) A‚f dt A ) k′(B′ - D′)

(24)

where B′ is the (upper triangular) matrix of elements {B′i,j} ) xixjR-2θ(xi/xj)∆xj and D′ is an N × N diagonal matrix with {D′i,i} ) xiR. Equation 24 describes N linear ODEs that can be easily and efficiently solved numerically, because the matrix A is triangular. An analytical solution15 in terms of eigenvalues/eigenvectors can also be given, after the unique diagonalization of A, although the calculation of a large number of eigenvalues, most


Figure 6. Simulation of the evolution of an initially monomodal distribution by the β-PDF method (solid line) and comparison with the method of classes (dashed line). Approximation with one β-PDF and three collocation points ([). Final time ) 3, k′) 1, R ) 2.5, binary β-daughter distribution θ(z), as in eq 23.

of them very close to each other, can be numerically difficult. Results are now presented for power-law and binary β-daughter distributions, for monomodal and bimodal initial distributions, and a variable number of β-PDFs in the approximation of eq 10. Note that the power law of eq 18 actually implements an inverse dependence on z because its exponent is always negative; this causes a large number of fine particles to result from breakage events, and the smaller the parameter γ, the smaller the daughters and the higher their number. Such a process induces an extensive broadening of the initial distribution. This is not the case for the binary β-distribution, which implements a breakage process that halves the particle volume, on average, with some variance about this mean value. Eventually, a second population of finer particles gradually emerges and shifts toward smaller volumes. It is then easier to introduce the method starting with the simulation of a breakage process that follow the binary β-daughter distribution. Figure 6 shows the evolution of an initially monomodal distribution f(x,0) ) f0(x) as simulated by two methods, the one proposed here (β-PDFs, with just one β-function in the approximation of eq 10) and the reference one, the method of classes applied to 200 sizes. The initial distribution was generated as β(x,15,5). Other parameter values are reported in the figure caption. Neither the units nor the value of the final time are relevant; both are related to the units and value of k′. A longer evolution of the distribution can be simulated by either increasing the final time or decreasing the kinetic constant k′. The final time for the presentation of the results is chosen to show the transient, most critical configurations. As expected, Figure 6 shows that a single β-PDF cannot reproduce the multimodal transients. Although the first moments for the instantaneous distributions calculated by the two methods, at the same time, remain fairly close (0.6% error on the first moment and 7.5% on the second moment at t ) 3), there is a structural deficiency of the single β-PDF that cannot be accommodated by a different arrangement or number of collocation points. By simply adding a second component

Figure 7. Same conditions of Figure 6, with two subdistributions in eq 10 and the minimum number of collocation points (five). Final time ) 5.

Figure 8. Same conditions of Figure 7, with one additional collocation point. Final time ) 10.

to the approximation of eq 10, we obtain the results of Figure 7. The minimum number of collocation points (five) uniformly distributed in the 0-1 range has been used. Minor local discrepancies can be identified. They disappear upon adjustment of the positions of the collocation points. An optimization criterion for the selection of collocation points should be developed. More easily, a single extra collocation point makes the distributions predicted by the two methods barely distinguishable, as shown in Figure 8, even when collocation points uniformly distributed in the whole range are used. Still, we are solving only six ODEs, and the maximum error throughout the evolution of the distributions on the evaluation of the first 15 moments with respect to the solution by classes (200 ODEs) is less than 4.6% for the 15th and 1.3% for the first. The single β-PDFs, β1 and β2, that combined give the result of Figure 8 are shown in Figure 9. One clearly sees that β1 describes almost alone the initial distribution and its initial evolution, slightly shifting toward smaller sizes. Then, β2 progressively emerges to dominate the approximation of the later stages of the total distribution. Circumstances in which one or the other

2656


Figure 9. Evolution of the two (weighted) components of g(x,t). Their sum results in the solid lines of Figure 8. Collocation points are also shown.

Figure 10. Simulation of the evolution of an initially monomodal distribution. Approximation with two β-PDFs. Final time ) 1, k′ ) 1, R ) 2, power-law daughter distribution θ(z), as in eq 18 with γ ) 0.1.

component is completely absent are never observed. At each stage of the evolution, both components are present, although one can be more relevant while the other acts as a correction. Interestingly, the subdistributions are accurately predicted even when the collocation points do not entirely cover the size domain where the distributions are located. The power-law daughter distribution behaves quite differently from the binary β-distribution, as explained above. It produces a larger amount of fines, particularly for small values of the parameter γ. The simplest case of uniform breakage (γ ) 1) results in an evolution comparable to that shown in Figures 6-8, and the same comments on the limitation of a single β-PDF and the requirement of a sufficient number or optimal choice of collocations points apply. Quite different results are obtained with small values of γ. The massive production of fines is reflected by a long tail toward the small sizes, as shown in Figure 10. Such a skewness of the distribution cannot be accommodated by a single β-PDF, even when the collocation points are increased in number or redistributed. A single β-PDF shifts to smaller x values to account for the dramatic decrease of the average size

Figure 11. Evolution of the two (weighted) components of g(x,t). Their sum results in the solid lines of Figure 10. Collocation points are also shown.

Figure 12. Simulation of the evolution of an initially bimomodal distribution. Approximation with three β-PDFs. Final time ) 30, k′ ) 1, R ) 2.5, binary β-daughter distribution θ(z), as in eq 23.

caused by such a breakage process, but it cannot fit such an asymmetric distribution. Even though the resulting distribution is not exactly multimodal, it clearly contains a developing population of fines that, together with the original population of larger particles, results in a very broad distribution. Such a process can be described by the suggested method with two subdistributions, where the second accounts for the increasing amount of fines. Results for this case are shown in the same Figure 10, where the role of the second subdistribution in describing the production of smaller particles is evident. It is of some interest to inspect the shape of the two components, which is shown in Figure 11. The second component, β2, reflects the flexibility of the β-PDF in representing distributions that are very skewed toward small sizes. As a final example, we apply the method to a higher multimodality. Assuming that the initial distribution is strictly bimodal, its evolution is expected to display a higher degree of multimodality. Figure 12 shows the evolution of an initial distribution f(x,0) ) f0(x) generated as the function 0.9β(x,15,5) + 0.1β(x,50,100), i.e., 90% is provided by the same initial distribution f0(x) used in the previous calculations, with the addition of


Figure 13. Evolution of the three (weighted) components of g(x,t). Their sum results in the solid lines of Figure 12. Collocation points are also shown.

10% of a population of smaller particles, centered at x ) 50/150 ) 1/3. The method was applied with three β-PDFs, as suggested by the considerations above and by the observation of the solution. Both of the initial distributions are affected by the breakage process, which moves mass from the region of larger particles in each distribution to the smallerparticle side. The intermediate distribution, however, also receives particles from the larger population component, broadening on both sides. Overall, the average size decreases, and the initial differences tend to vanish in a single, fairly broad monomodal distribution. During transients, broad asymmetrical bimodal distributions are observed. The evolution of the three components β1, β2, and β3 is shown in Figure 13. The growth of the final distribution of fines is described by the third component, while the first two β-PDFs, describing the initial distribution, progressively disappear, partially shifting to smaller sizes. Conclusions In this work, we formulated a method to predict the evolution of generic populations undergoing pure breakage. Prediction is achieved by solving the PBE for given initial conditions and breakage functions. The method is specifically conceived for distributions that can show multimodality of any order. It can also be used for pure monomodal distributions, but we do not see any specific advantage in this case with respect to other known techniques. We believe that the method suggested takes advantage of some original ideas. One is the suggestion of tracking the evolution of the particle size distribution through an analytical PDF with a suitable shape, adjusting its parameters. This provides the opportunity to have a complete description of the distribution at any time, instead of specifying only limited values of it (classes) or some features (moments). In addition, a careful selection of the approximating PDF can lead to simplifications in the evaluation of the birth and death functions in the PBE, allowing explicit evaluation of the integral. At the same time, the use of analytical PDFs reduces the amount of information required to describe a complex function, such as a distribution, to the number of parameters in the analytical expression of

the PDF(s). Even though methods based on standard or problem-specific PDFs have already been presented,5,7 we introduced here two different concepts: the simplest is using β-PDFs as approximating functions, which are very flexible and analytically amenable to many simplifications in the evaluation of the integrals, with several breakage functions. The second original concept is using many PDFs for the approximation, with unknown functions as their parameters, instead of developing polynomial expansions based on a single PDF, as in refs 5 and 7. This latter concept is crucial to providing a physical basis to the description and simulation of multimodal distributions, because it suggests that multimodality can be the result of the coexistence of several subpopulations in the sample. The concepts above have been developed both (1) demonstrating that the parameters in the standard PDF can be used directly to fit the experimental data and (2) deriving the equations for the variation of the PDF parameters from the original PBE, with extensive analytical pretreatment. The results presented allow one to conclude that the approach can be very efficient. Most of the results are amenable to local improvement using additional PDFs in the weighted approximant of eq 10, or using a more appropriate distribution of collocation points along the size interval, or even formulating a different solution scheme when the PBE has been translated into a single ODE, such as eq 22. We also believe that the method can be extended to include the case of size enlargement processes whenever an upper bound to the distribution can be identified, thereby allowing the size to be scaled to a 0-1 range. Acknowledgment This work was partially supported by the FIRB contract N. RBAU01NZH7 and the PRIN contract N. 2003093941. Literature Cited (1) Ray, W. H. On the mathematical modeling of polymerization reactors. J. Macromol. Sci.-Rev. Macromol. Chem. 1972, C8, 1. (2) Ramkrishna, D. Population Balances: Theory and Application to Particulate Systems in Engineering; Academic Press: San Diego, 2000. (3) Metha, K.; Madras, G. Dynamics of molecular weight distributions for polymer scission. AIChE J. 2001, 47, 2539. (4) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes: Analysis and Techniques of Continuous Crystallization; Academic Press: New York, 1988. (5) Hamilton, R. A.; Curtis, J. S.; Ramkrishna, D. Beyond lognormal distributions: Hermite spectra for solving population balances. AIChE J. 2003, 49 (9), 2328. (6) Kostoglous, M.; Karabelas, A. J. On the breakage problem with a homogeneous erosion type kernel. J. Phys. A: Math. Gen. 2001, 34, 1725. (7) Canu, P.; Ray, W. H. Discrete weighted residual methods applied to polymerization reactions. Comput. Chem. Eng. 1991, 15 (8), 549. (8) Garcia, F.; Le Bolay, N.; Frances, C. Changes of surface and volume properties of calcite during a batch wet grinding process. Chem. Eng. J. 2002, 85, 177. (9) Bel Fadhel, H.; Frances, C.; Mamourian, A. Investigations on ultra-fine grinding of titanium dioxide in a stirred media mill. Powder Technol. 1999, 105, 362. (10) Peleg, M.; Normand, M. D. Simulation of size reduction and enlargement processes by a modified version of the Beta distribution function. AIChE J. 1986, 32 (11), 1928. (11) Popplewell, L. M.; Campanella, O. H.; Peleg, M. Simulation of bimodal size distribution in aggregation and disintegration processes. Chem. Eng. Prog. 1989, 85, 56.

2658


(12) Canu, P.; Strumendo, M. Efficient treatment of PBE for size reduction processes. In Proceedings of PARTEC 2004; Pratsinis, S. E., Schulz, H., Strobel, R., Schreglmann, C., Eds.; Nu¨rnbergMesse GmbH: Nu¨rnberg, Germany, 2004; p 224. (13) Molina-Boisseau, S.; Le Bolay, N. Fine grinding of polymers in a vibrated bead mill.Powder Technol. 1999, 105, 321. (14) Bel Fadhel, H.; Frances, C. Wet batch grinding of alumina hydrate in a stirred bead mill. Powder Technol. 2001, 119, 257. (15) Berthiaux, H.; Dodds, J. A new estimation technique for the determination of breakage and selection parameters in batch grinding. Powder Technol. 1997, 94, 173. (16) Lensu, M. Correlation between fragment sizes in sequential fragmentation. J. Phys. A: Math. Gen. 1997, 30, 7501. (17) Kostoglous, M.; Dovas, S.; Karabelas, A. J. On the steadystate size distributions of dispersions in breakage processes. Chem. Eng. Sci. 1997, 52, 1285. (18) Kapur, P. C. Self-preserving size spectra of comminuted particles. Chem. Eng. Sci. 1972, 27, 425. (19) Austin, L. G.; Luckie, P. T. Estimation of nonnormalized breakage distribution parameters from batch grinding. Powder Technol. 1972, 5, 267. (20) Jones, W. P.; Whitelaw, J. H. Calculation Methods for Reacting Turbulent Flows: A Review. Combust. Flame 1982, 48, 1. (21) Vanni, M. Discretization procedure for the breakage equation. AIChE J. 1999, 45, 916. (22) Diemer, R. B.; Olson J. H. A moment methodology for coagulation and breakage problems: Part 2sGeneralized daughter distribution functions. Chem. Eng. Sci. 2002, 57, 4187.

(23) McCoy, B. J.; Madras, G. Analytical solution for a population balance equation with aggregation and fragmentation. Chem. Eng. Sci, 2003, 58, 3049. (24) Elhanbaly, A. On the solution of the integro-differential fragmentation equation with continuous mass loss. J. Phys. A: Math. Gen. 2003, 36, 8311. (25) Teymour, F.; Campbell, J. D. Analysis of the dynamics of gelation in polymerization reactors using the “numerical fractionation” technique. Macromolecules 1994, 27, 2460. (26) Hulburt, H. M.; Katz, S. Some problems in particle technology. A statistical mechanical formulation. Chem. Eng. Sci. 1964, 19, 555. (27) Pladis, P.; Kiparissides, C. A Comprehensive Model for the Calculation of Molecular Weight-Long-Chain Branching Distribution in Free-Radical Polymerizations. Chem. Eng. Sci. 1998, 53, 3315. (28) Kruse, T. M.; Wong, H.-W.; Broadbelt, L. J. Modeling the Evolution of the Full Polystyrene Molecular Weight Distribution during Polystyrene Pyrolysis. Ind. Eng. Chem. Res. 2003, 42, 2722. (29) Hill, P. J.; Ng, K. M. New Discretization Procedure for the Breakage Equation. AIChE J. 1995, 41 (5), 1204.

Received for review April 29, 2004 Revised manuscript received September 4, 2004 Accepted September 29, 2004 IE049652C

Prediction of Multimodal Distributions in Breakage Processes

Recommend Documents