4380
Ind. Eng. Chem. Res. 2004, 43, 4380-4393
Integrated Framework for the Numerical Solution of Multicomponent Population Balances. 1. Formulation, Representation, and Growth Mechanisms Darren D. Obrigkeit and Gregory J. McRae* Department of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Room 66-372, Cambridge, Massachusetts 02139
Current numerical solutions of population balance problems including growth are limited by the Courant condition time step constraint. As a result, particle size ranges must be restricted and limited growth rate mechanisms can be solved. A thorough analysis of growth rate mechanisms and dynamics is presented that relates the growth rate mechanisms and particle size ranges to the time step limit. These relationships reveal a method for producing the optimal time step for any growth rate model, which can increase the time step by a factor of 107. When these methods are combined with previous approaches, new hybrid methods are developed that allow the user to flexibly choose the optimal time step given the resolution requirements of the problem, producing solutions that are 104-107 times faster than conventional solutions. The numerical accuracy of these solutions compares favorably against analytical solutions, demonstrating their applicability across a wide range of problems. Introduction Dispersed-phase systems are the cornerstone for many classes of physicochemical processes, including aerosols, crystallization, combustion, and precipitation. These processes have critical end applications ranging from air pollution1 to polymers,2 protein characterization,3 and ceramic biomaterials.4 Despite the widespread importance of these processes, their dynamics are poorly understood. The complex nature of competing coagulation, nucleation, and growth mechanisms over many orders of magnitude in particle size results in complex dynamics that vary widely in time and length scales. Coupled with the limitations of particle measurement systems, it is virtually impossible to resolve the interactions of these mechanisms through experimentation. A number of barriers likewise hinder the ability of models to describe population balance system dynamics. The combination of multicomponent systems, particle size ranges covering up to 12 or more orders of magnitude, and the partial integrodifferential equation systems necessary to describe coagulation and agglomeration produces essentially intractable problems. These barriers limit the types of mechanisms and particle size ranges that can be described by current numerical solutions. The highly restricted set of cases where analytical solutions are applicable to population balance systems precludes their use as a general modeling tool; the only true utility of analytical solutions5 is to confirm the accuracy of numerical solutions. This work introduces new advances that for the first time enable dynamic simulation of complex multicomponent systems over a large particle size range. In particular, new representations have been developed for multicomponent systems that are compact yet still capture variations in particle composition, allowing solutions that are up to 5 orders of magnitude faster than conventional * To whom correspondence should be addressed. Tel.: (617) 253-6564. Fax: (617) 258-0546. E-mail:
[email protected].
methods. New transformations have been developed to convert data from traditional and newly developed representations; a general formulation has been developed to derive these transformations between any arbitrary set of representations. A new analysis of growth mechanisms has been conducted, revealing new insights into the dynamics of particle growth and leading to the discovery of new numerical solution methods covering the full range of particle sizes necessary to accurately reproduce particle size distributions. These methods produce solution times up to 107 times faster than conventional methods while allowing the user the flexibility to choose the solution resolution over the particle size range. This work reveals these new advances in transformations, growth mechanism analysis, and new numerical methods that enable a description of the entire particle size range; methods for formulating multicomponent representations, their transformations, and a detailed description of their implementation in a dynamic numerical solution will follow in part 2 of this series of papers. The inability of current numerical solutions to cover the entire particle size range (see Figure 1) often leads to large errors in the number density and size of particles as well as in the overall mass balance for the system. This shortfall represents a critical barrier in developing accurate models for a wide range of particulate processes; a survey of recent solution methods reveals that only in the special case of a volume reaction-limited growth rate mechanism can more than 2 orders of magnitude in particle size be modeled6-11 (see Table 1). Attempts have been made to circumvent this issue by using moments to describe the shape of the number density distribution. The simplest moment method approach assumes a monodisperse distribution,12 while more sophisticated approaches assume a predetermined particle size distribution function, such as log-normal.13 More elaborate approaches describe multimodal particle size distributions by dividing the overall distribution into a number of individual mono-
10.1021/ie020548+ CCC: $27.50 © 2004 American Chemical Society Published on Web 06/04/2004
Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4381
distribution resolution and can be solved up to 107 times faster than conventional methods are introduced. These methods are verified against analytical results and are shown to be in good agreement. Representation of Population Balance Systems
Figure 1. Comparison of typical particle size ranges in a number of population balance systems with particle size range described by current numerical solution methods. Note that, in general, these solution methods are restricted to a small portion of the total particle size range.
disperse14 or log-normal modes.15 Whitby et al. have developed a general framework for multimodal problems that manages the creation, elimination, and consolidation of modes.16 Other methods have sought to relax these assumptions restricting particle size distribution form by using expansions of orthogonal basis functions. Typically, the basis set is constructed from Laguerre polynomials because their characteristic domain [0, ∞] naturally covers the entire possible particle size range.17,18 Nevertheless, these methods are inefficient and produce poor representations of the particle size distribution even when 20-50 Laguerre expansion terms are used.19-21 Still others have employed closure methods to estimate the nonpositive and noninteger moments needed to produce a closed set of equations for many applications.20,22 While these methods are very efficient and do not restrict the form of the particle size distribution, they do not provide sufficient information to uniquely reconstruct the particle size distribution and are therefore only acceptable in applications where distribution moments are the key performance metric. Diemer and Olson21 propose basis functions combining steady-state solutions in the small particle size range with generalized exponential basis functions to reconstruct particle size distributions and demonstrate close reproduction of coagulation-breakage problems, while Wulkow19 uses an adaptively placed variable-order local Chebyshev polynomial expansion to construct a particle size distribution. Even though these methods allow full reconstruction of the particle size distribution, they do not address the issue of varying growth rate mechanisms. This is a critical barrier to the accurate implementation of population balance models and their ability to describe growth mechanisms in the design and optimization of particulate processes. This work presents a framework for solving population balance systems with arbitrary growth rate mechanisms over large particle size ranges. First, a brief background on population balance systems and number density representation is presented, followed by a more detailed discourse on growth rate mechanisms, their numerical stability requirements, and the resulting dynamics of numerical solution methods. New scaling methods that allow the flexible selection of particle size
The number density n(φ,t) is the fundamental function used to describe the number of particles at each size in an ensemble of particles; it is defined such that n(φ,t) dφ is the number of particles occupying some amount of space in the size range [φ, φ + dφ] at time t (Figure 2). As a result, the units of n(φ,t) dφ are #/m3. All density functions express a frequency of occurrence over some range of independent variable values. For the special case of the number density, the density function represents the frequency of occurrence of particles, while a probability density function23,24 represents the frequency of occurrence of events. Table 2 summarizes the similarities between number density functions and probability density functions. Because n(φ) represents the incremental increase in the total number of particles N as the independent variable is increased in the cumulative number density integral, n(φ) is also commonly expressed as dN/dφ. Also, note that φ can represent any property of particles in a number density distribution, including mass, volume, crystal length, or composition. As φ takes on different units, the units of n(φ) will also change such that n(φ) dφ will always have units of number of particles/unit volume. Table 3 shows several examples of the various units that φ can take on and the corresponding units for n(φ) as well as the dN/dφ expression. It is often necessary to derive transformations that translate data from one form of n(φ) to another. For instance, experimental systems25 typically measure particle diameter (Dp) while numerical solutions often use mass, m, as the independent distribution variable. To compare these data, a transformation is needed from n(Dp) to n(m). The key governing relation in this transformation is that corresponding segments in each number density representation must contain the same number of particles:
∫mm n(m) dm ) ∫DD 2
p2
1
p1
n(Dp) dDp
(1)
where Dp1 ) f(m1) and Dp2 ) f(m2) are related by the formula
Dp(m) )
3
FπDp 6m 1/3 S m(Dp) ) Fπ 6
( )
(2)
Equation 1 can also be written in differential form:
n(m) dm ) n(Dp) dDp
(3)
which in combination with eq 2 yields the correct transformation
n(m) ) n(Dp)
dDp dm
dDp 1 6 -2/3 ) m dm 3 Fπ
( )
2 ∴n(m) ) n(Dp) FπDp2
(4)
4382 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 Table 1. Comparison of Particle Size Ranges and Growth Rate Mechanisms for Conventional Numerical Solutionsa system
growth law
particle size range (orders of magnitude)
ref
crystallizer crystallizer crystallizer aerosol polymerization polymerization aerosol aerosol
diffusion constant constant surface reaction volume reaction constant volume reaction volume reaction
1 2 2 2 1 5 5 9
Pantelides and Oh, 1996 Chung et al., 1999 Litster et al., 1995 Katoshevski and Seinfeld, 1997 Alvarez et al., 1992 Wulkow, 1996 Kim and Seinfeld, 1990 Pilinis, 1990
a In general, solutions are incapable of describing more than 2 orders of magnitude in particle size except in the special case of a volume-limited reaction growth mechanism.
Table 2. Similarities between Probability Density Functions and Number Density Functions as Well as Integrals over Respective Intervals of These Functions number density function
probability function
interpretation
∫ n(O) dO
total particles in the interval (a, b)
∫
b
n(O) dO
cumulative number of particles of size
