Solution of Population Balance Equations by the Finite Size Domain

5222

Ind. Eng. Chem. Res. 2010, 49, 5222–5230

Solution of Population Balance Equations by the Finite Size Domain Complete Set of Trial Functions Method of Moments (FCMOM) for Inhomogeneous Systems Matteo Strumendo* and Hamid Arastoopour Department of Chemical and Biological Engineering, Illinois Institute of Technology, Chicago, Illinois 60616

The FCMOM (finite size domain complete set of trial functions method of moments) is an efficient and accurate numerical technique to solve monovariate and bivariate population balance equations. It was previously formulated for homogeneous systems. In this paper, the FCMOM approach is extended to solve monovariate population balance equations for inhomogeneous (spatially not uniform) systems. In the FCMOM, the method of moments is formulated in a finite domain of the internal coordinates and the particle size distribution function is represented as a truncated series expansion by a complete system of orthonormal functions. The FCMOM is extended to inhomogeneous systems assuming that the particle-phase convective velocity is independent of the internal variables (particle size). The method is illustrated by applications to particle diffusion and to particle convection. In the case of particle convection, a gas-solid dilute flow in a pipe was simulated and the FCMOM equations were coupled with the governing equations (mass and momentum balances) of the gas phase. 1. Introduction Population balance equations (PBE) are widely used in chemical engineering to model processes involving one or more particulate phases. Applications of PBE range from fluidized beds to droplet dispersions, from cell population evolution to crystallization, from size reduction to particle aggregation processes, from polymerization reactions to aerosol formation.1 In the area of fluidization of gas-solid systems, in several relevant processes such as catalytic cracking, coal combustion, and coal gasification the particle polydispersity affects both the gas/solid fluid dynamics and the heterogeneous reaction rates; conversely, the particle size distribution is determined, among other factors, by the reaction rates and by the gas/solid fluid dynamics (for example, in a coal gasifier one of the particle fragmentation mechanisms is through thermal stresses). The book of Ramkrishna1 provides a comprehensive review of the subject of PBE in terms of formulation of PBE, solution techniques, PBE applications, and theoretical considerations. In recent years, considerable attention has been paid to a specific subject in the area of PBE, i.e., the development of numerical solution techniques based on the method of moments (MOM) which demand a small computational effort. Such a requirement is important in the perspective of coupling the PBE with the transport equations (mass, momentum, energy, and species) of the fluid and particulate phases. More specifically, three variations of the MOM, capable of providing numerical solutions of PBE using efficient algorithms, were developed, namely, the quadrature method of moments (QMOM), the direct quadrature method of moments (DQMOM), and the finite size domain complete set of trial functions method of moments (FCMOM). The QMOM was proposed by McGraw2 and then was validated (for monovariate PBE) by Barrett and Webb3 and by Marchisio et al.45 In the QMOM, no explicit assumption is made regarding the shape of the size distribution function and the integrals appearing in the moments equations are computed numerically by means of quadrature formulas. * To whom correspondence should be addressed. Tel.: 001 312 567 3066. Fax: 001 312 567 8874. E-mail: [email protected].

The DQMOM, proposed by Marchisio and Fox,6 differs from the QMOM because the abscissas and the weights of the quadrature approximation are tracked directly (rather than the moments, as in the QMOM). Furthermore, an explicit expression for the particle distribution function is given in terms of a summation of Dirac δ functions. Additionally, using the DQMOM, the dependency of the particle velocities with respect to the internal variables is accounted for. The FCMOM was presented and validated for monovariate PBE in homogeneous conditions by Strumendo and Arastoopour7 and extended to bivariate PBE by Strumendo and Arastoopour.8 Some advantages of the FCMOM (for monovariate and bivariate PBE) are the following: (1) it provides accurate values both for the moments and for the reconstructed particle size distribution (or reconstructed bivariate particle distribution function); (2) it is computationally efficient (even in bivariate applications); (3) the internal variables always assume finite values, and the domains of the internal variables are clearly defined within boundaries, thus assuring that the values assumed by the internal variables are always realistic. For a more detailed discussion regarding the different variations of the MOM and their comparison, see Strumendo and Arastoopour.7,8 The goal of this work is to validate the FCMOM in the solution of monovariate PBE for inhomogeneous systems. The general structure of a monovariate PBE includes the following terms (Marchisio and Fox6 and Hulburt and Katz9): (1) accumulation rate; (2) terms accounting for spatially inhomogeneous conditions; (3) source terms. The terms accounting for spatially inhomogeneous conditions include convection and spatial diffusion terms. The source terms include terms due to nucleation, particle growth and/or dissolution, particle aggregation and breakage. In this paper, the focus is on the treatment (in the FCMOM framework) of the convection and spatial diffusion terms; instead, the source terms are neglected, because their treatment was thoroughly investigated by Strumendo and Arastoopour.7 When dealing with PBE for inhomogeneous systems, the convection term can be modeled making the assumption that the particle-phase velocity is only a function of the external

10.1021/ie901407x  2010 American Chemical Society Published on Web 04/21/2010

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

coordinates (time and spatial coordinates) or, conversely, that the particle-phase velocity depends both on the external coordinates and on the internal coordinate (size). The former assumptionsi.e., particle-phase velocity independent of the particle sizesis justified when considering light particles in a dilute flow, while it may be not acceptable for massive particles. Presently, a method capable of taking into account the dependency of the particle velocities with respect to the internal variables is the DQMOM (Marchisio and Fox6); this capability of the DQMOM was demonstrated in simulations of polydisperse gas-solid systems in dense fluidized beds (Fan et al.10 and Fan and Fox11). In this work, the assumption that the particle-phase velocity is only a function of the external coordinates is still retained; the extension of the FCMOM to include particle-phase velocities dependent also on the particle size is currently under investigation by the authors. The paper is organized as follows. In the second section, the governing equations of the FCMOM technique for inhomogeneous processes are presented and discussed. In the third and fourth sections, applications to spatial diffusion and convection are considered separately. In the third section, the FCMOM is applied to the diffusion equation. In the fourth section, particle convection in a gas-solid pipe flow is considered; the FCMOM equations are coupled in a bidimensional geometry with the equations governing the gas-phase fluid dynamics. In the fourth section, the particle-phase velocity is further assumed to be equal to the gas velocity. However, this is not a restriction of the methodology presented in this paper; the only restriction is, as already mentioned, that the particle-phase velocity is only a function of the external coordinates (not a function of the particle size). Assuming that the particle-phase velocity is equal to the gas velocity helps us to focus our attention on the FCMOM validation and to check carefully that the particle size distribution function is correctly computed.

∂Vp, j(t, x) f(ξ, t, x) ∂f(ξ, t, x) + ∂t ∂xj ∂ ∂f(ξ, t, x) D (ξ, t, x) ) 0 (2) ∂xj pt ∂xj

[

]

or equivalently

[

]

∂Dpt(ξ, t, x) ∂Vp, j(t, x) ∂f(ξ, t, x) + Vp, j(t, x) × + f(ξ, t, x) ∂t ∂xj ∂xj ∂f(ξ, t, x) ∂ ∂f(ξ, t, x) - Dpt(ξ, t, x) ) 0 (3) ∂xj ∂xj ∂xj

[

]

Using the following coordinate transformation ξ)

{ξ - [ξmin(t, x) + ξmax(t, x)]/2} [ξmax(t, x) - ξmin(t, x)]/2

(4)

the dimensionless internal variable ξj ranges in the domain [-1, +1]. The particle size distribution function can be defined in dimensionless terms as f′(ξ, t, x) ) f(ξ, t, x)/fsc where fsc is an appropriate scale factor. The dimensionless particle size distribution function f′ is expressed as a series expansion truncated after the first M terms: M-1

f′(ξ,t, x) =

∑ c (t, x) φ (ξ) n

(5)

n

n)0

In the series expansion 5, φn(ξj) are the orthonormal functions associated with the Legendre polynomials (Courant and Hilbert12) Pn(ξj) by φn(ξ) )

2. Extension of the FCMOM to Inhomogeneous Systems The FCMOM was presented and validated for monovariate PBE in homogeneous conditions by Strumendo and Arastoopour.7 In this section, such an approach is extended to monovariate PBE for inhomogeneous systems. The particles are characterized by one internal variable, i.e., the particle radius ξ, whose values range between a minimum value ξmin and a maximum value ξmax. The external coordinates are the time t and the spatial Cartesian coordinates x (or x, y, z). The particle size distribution is denoted by f(ξ,t,x). A general formulation of monovariate PBE is (McGraw,2 Marchisio and Fox,6 and Hulburt and Katz9)

5223


2n 2+ 1 P (ξ)

(6)

n

The first few Legendre polynomials are 3 1 P2(ξ) ) ξ2 - , 2 2 5 3 35 4 15 2 3 P3(ξ) ) ξ3 - ξ, P4(ξ) ) ξ ξ + 2 2 8 4 8

P0(ξ) ) 1,

P1(ξ) ) ξ,

In the series expansion 5, the coefficients of the series cn are given by cn )

∫

+1

-1

f′φn dξ

(7)

f′(ξ)i dξ

(8)

The moments are defined as ∂Vp, j(t, x) f(ξ, t, x) ∂f(ξ, t, x) + ∂t ∂xj ∂ ∂f(ξ, t, x) ∂G(ξ, t, x) f(ξ, t, x) + D (ξ, t, x) )∂xj pt ∂xj ∂ξ B(ξ, t, x) - D(ξ, t, x) (1)

[

]

On the left-hand side of eq 1, vp(t,x) is the particle-phase velocity and Dpt is the turbulent diffusivity. On the right-hand side of eq 1, G is the particle growth rate and B and D are the birth and death terms due to nucleation, breakage, and aggregation, etc. The source terms on the right-hand side of eq 1 due to particle growth, nucleation, breakage, and aggregation are neglected in this work (see Strumendo and Arastoopour7 for the treatment of such terms). Therefore, the PBE can be rewritten as

µi )

∫

+1

-1

and the coefficients cn can be expressed in terms of the moments by cn )




n

2n + 1 1 (2ν)! × (-1)n-ν 2 2n ν)0 [(2ν - n)!] 1 µ (9) [(n - ν)!][(ν)!] 2ν-n

∑

{

}

where the terms with negative moments order are to be omitted. The particle size distribution function f′ is therefore completely determined by the moments. The moments evolution equations are derived from eq 3, after applying the coordinate transformation (4):

5224

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

[∫

]

∂µi ∂Vp, j ∂µi +1 ∂f ′ i ∂ + µi + Vp, j D′ (ξ) dξ ) ∂t ∂xj ∂xj ∂xj -1 pt ∂xj -(MB + MBconv + MBdiff1 + MBdiff2 + MBdiff3) (10) where

(

)

1 ∂ξmin + ξmax ∆ξ ∂t 1 ∂∆ξ {[f+1 ′ - (-1)i+1f-1 ′ ] - (i + 1)µi} ∆ξ ∂t

MB ) -{[f+1 ′ - (-1)if-1 ′ ] - iµi-1}

(

( )

)

Vp, j ∂ξmax + ξmin ∆ξ ∂xj Vp, j ∂∆ξ ′ ] - (i + 1) µi} {[f1′ - (-1)i+1 f-1 ∆ξ ∂xj

MBconv ) - {[f1′ - (-1)if-1 ′ ] - iµi-1}

MBdiff1 )

MBdiff2 )

)

∫

(

2 ∂(ξmax + ξmin) ∆ξ ∂xj

∫

1 ∂ξmax + ξmin ∆ξ ∂xj

( )

∂Dpt ′ ∂f′ i (ξ) dξ + -1 ∂xj ∂ξ ′ ∂f′ i+1 +1 ∂Dpt 1 ∂∆ξ (ξ) dξ -1 ∆ξ ∂xj ∂xj ∂ξ +1

( )∫

1

-1

Dpt ′

∂2f ′ (ξ)i dξ + ∂xj ∂ξ

2 ∂∆ξ 1 ∂2f ′ D (ξ)i+1 dξ ′ ∆ξ ∂xj -1 pt ∂xj ∂ξ 1 2 ∂(ξmax + ξmin) 2 1 ∂2f ′ i D (ξ) dξ ′ -1 pt ∆ξ ∂xj ∂ξ2 1 2 ∂(ξmax + ξmin) ∂∆ξ 1 ∂2f′ i+1 2 D (ξ) dξ ′ pt ∆ξ ∂xj ∂xj -1 ∂ξ2 1 2 ∂∆ξ 2 1 ∂2f ′ Dpt ′ 2 (ξ)i+2 dξ -1 ∆ξ ∂xj ∂ξ 2 ∂∆ξ ∂(ξmax + ξmin) 1 ∂f ′ i D′ (ξ) dξ + 2 ∂x -1 pt ∂x ∂ξ (∆ξ) j j 2 1 ∂ (ξmax + ξmin) 1 ∂f ′ i D′ (ξ) dξ 2 -1 pt ∆ξ ∂ξ ∂x

( )[ ( ) ( )(

[

∫

]∫

∫

)∫

j 2

∫

∫

( ) ]∫

∂ ∆ξ 2 ∂∆ξ 1 + 2 ∆ξ ∆ξ ∂xj ∂xj

(

)∫

2

∂f ′ i+1 D′ (ξ) dξ -1 pt ∂ξ

[ ( ( )(

)] ∫ ) ∫ (

up, j,min(t, x) ) up, j(ξ ) ξmin, t, x) ) Dpt(ξ ) ξmin, t, x) ∂f Vp, j(t, x) (ξ ) ξmin, t, x) f(ξ ) ξmin, t, x) ∂xj

[

up, j,max(t, x) ) up, j(ξ ) ξmax, t, x) ) Dpt(ξ ) ξmax, t, x) ∂f Vp, j(t, x) (ξ ) ξmax, t, x) f(ξ ) ξmax, t, x) ∂xj

[

)∫

′ , jf-1 ′ denote the values of jf ′ at ξj ) +1 and ξj ) In eq 10, jf+1 ′ (ξj,t,x) ) Dpt(ξ,t,x). -1; ∆ξ is equal to (ξmax - ξmin) and Dpt The structure of eq 10 is as follows. On the left-hand side, there are terms which are present even when ξmin and ξmax are constant and homogeneous; specifically, the first term is the moments accumulation rate, the second and third terms represent convection, and the fourth term represents particle diffusion. The second term is zero if the particle phase is considered incompressible.

]

]

The moving boundary condition for ξmin is then given by ∂ξmin ) ∂t

∑

i)x,y,z

(

-

)

∂ξmin u ∂xi p,i,min

and the moving boundary condition for ξmax is

1


1 ∂Dpt ∂f ′ 1 ∂ξmax + ξmin MBdiff3 ) (ξ)i dξ + -1 ∆ξ ∂xj ∂ξ ∂xj
1 ∂∆ξ 1 ∂Dpt ∂f ′ i+1 (ξ) dξ ∆ξ ∂xj -1 ∂ξ ∂xj
1 ∂ξmax + ξmin 2 1 ∂f ′ ∂Dpt i (ξ) dξ -1 ∆ξ ∂xj ∂ξ ∂ξ
1 2 ∂ξmax + ξmin ∂∆ξ 1 ∂f ′ ∂Dpt i+1 (ξ) dξ 2 ∆ξ ∂xj ∂xj -1 ∂ξ ∂ξ
1 ∂∆ξ 2 1 ∂f ′ ∂Dpt i+2 (ξ) dξ -1 ∆ξ ∂xj ∂ξ ∂ξ

∫

On the right-hand side of eq 10, there are the terms MB, MBconv, MBdiff1, MBdiff2, and MBdiff3 which are due to the coordinate transformation (4): MB accounts for the temporal variation of the boundaries and is present also in homogeneous processes (Strumendo and Arastoopour7). MBconv accounts for the spatial variation of the boundaries in the presence of convection. MBdiff1, MBdiff2, and MBdiff3 account for the spatial variation of the boundaries in the presence of particle diffusion. MBdiff1 is zero if the turbulent diffusivity is homogeneous, and MBdiff3 is zero if the turbulent diffusivity is assumed to be a function only of the external coordinates t, x. Moments evolution equations must be coupled with the moving boundary conditions providing the governing equations for ξmin(t,x) and ξmax(t,x). When the source terms (due to particle growth, nucleation, aggregation, and breakage) in eq 1 are neglected, ξmin(t,x) and ξmax(t,x) change because of the convective or diffusive particle movement in the presence of spatially inhomogeneous conditions. The definition of the moving boundary conditions for ξmin(t,x) and ξmax(t,x) is based on the evaluation of the following: (a) the spatial derivatives of ξmin(t,x) and ξmax(t,x); (b) the velocities up,min(t,x) and up,max(t,x) of the particles whose sizes are equal to ξmin(t,x) and to ξmax(t,x), respectively. The velocities up,min and up,max are expressed as the sum of the convective and of the diffusive contributions:

∂ξmax ) ∂t

∑

i)x,y,z

(

-

∂ξmax u ∂xi p,i,max

)

The moments evolution eq 10 coupled with the moving boundaries conditions form a set of partial differential equations whose variables are the moments µi(t,x), ξmin(t,x), and ξmax(t,x). The initial conditions for the moments are computed from the initial particle size distribution function by eq 8. 3. Applications to Particle Diffusion In this section the proposed approach is applied to pure diffusion processes. In the first application (subsection 3.1), the particle turbulent diffusivity is assumed to be size-independent. In subsection 3.2, the particle turbulent diffusivity is assumed to be size-dependent. According to the Tchen’s theory (Baldyga and Orciuch13 and Hinze14), the particle turbulent diffusivity is size-independent and equal to the fluid turbulent diffusivity for long diffusion times (compared to the Lagrangian integral time scale for the particles); instead it is size-dependent for short and intermediate diffusion times (compared to the Lagrangian integral time scale for the particles). In subsection 3.2 the limiting case of the Tchen’s theory for short diffusion times is considered.

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

5225

In the application presented in subsection 3.2 the moments equations are not closed, and therefore the FCMOM closure (eqs 5 and 9) is used to solve the moments equations. The capability of the FCMOM of predicting accurately the particle size distribution evolution for spatially homogeneous systems when the moments equations are not closed was demonstrated by Strumendo and Arastoopour;7 the same capability is here discussed for inhomogeneous systems. 3.1. Size-Independent Particle Diffusivity. In this subsection the particle turbulent diffusivity is assumed to be sizeindependent. Further, the particle turbulent diffusivity is assumed to be (spatially) homogeneous and constant (with respect to the time). In such conditions, and if a one-dimensional Cartesian geometry is used, the PBE (2) becomes the diffusion equation: ∂f ∂2f - Dpt 2 ) 0 ∂t ∂x

(11) Figure 1. PSD vs particle radius: initial PSD at different locations for the pure particle diffusion problem.

Correspondingly, in the evolution equations of the moments (10) all the convective terms and all the terms with the spatial derivative of the diffusivity are neglected. The domain of the Cartesian coordinate x is between 0 and L. The problem is solved with Dirichlet boundary conditions; i.e., the particle size distribution is prescribed for each time at x ) 0 and at x ) L. It is assumed that the initial particle size distribution (PSD) is given by the following distribution function: f(ξ, 0, x) ) a

{

}

q(x) [ξ - ξmin(0, x)] × [ξmax(0, x) - ξmin(0, x)] [ξmax(0, x) - ξ] [ξmax(0, x) - ξmin(0, x)]

{

}

s(x)

(12)

where a is a dimensional constant and the functions q(x) and s(x) are expressed as linear functions of x (q1, q2, s1, and s2 are constant values): q(x) )

x q + q2, L 1

x s(x) ) - s1 + s2 L

In the simulations, the following conditions and values were used. The constant values in eq 12 are q1 ) 6, q2 ) 2, s1 ) 6, and s2 ) 8. The boundary conditions are f(ξ,t,0) ) f(ξ,0,0) and f(ξ,t,1) ) f(ξ,0,1); i.e., at the boundaries the particle size distribution is kept equal to the initial particle size distribution function. Finally, at the initial time, ξmin(0,x) and ξmax(0,x) are assumed to be homogeneous (constant with respect to x) and equal to 10 and 70 µm, respectively. The simulation results are shown in Figures 1-3. Figure 1 shows the initial PSD function at x/L ) 0, 0.25, 0.5, 0.75, 1. Figures 2 and 3 show the numerical results at the final time in two different locations, compared with the exact solution (the exact solution is obtained using the Fourier series method). In Figure 2 the comparison between the numerical solution obtained by the FCMOM and the exact solution at x/L ) 0.25 is plotted; Figure 3 shows the same plot for x/L ) 0.5. In both locations the comparison between the numerical solution, obtained with 11 moments (M ) 11), and the exact solution is excellent. Further, the computational effort required to perform the simulations is low (less than 1 s). 3.2. Size-Dependent Particle Diffusivity. The one-dimensional diffusion process discussed in subsection 3.1 is considered in this subsection as well. However, here a size-dependent particle diffusivity is considered.

Figure 2. PSD vs particle radius: comparison between the numerical solution by the FCMOM and the exact solution for pure particle diffusion (sizeindependent diffusivity). The location is at x/L ) 0.25 and the final time is L2/Dpt. The numerical solution is obtained with M ) 11.

More specifically, the particle diffusivity is assumed to follow the Tchen’s expression for short diffusion times (Baldyga and Orciuch13 and Hinze14): 1 + b2(τp /TL) Dpt(ξ) ) Dft 1 + (τp /TL)

(13)

In eq 13, Dft is the fluid turbulent diffusivity and the Lagrangian integral time scale for the fluid motion TL can be computed from the turbulent kinetic energy kt and the turbulent kinetic energy dissipation rate εt using the expression: TL = 0.8

kt 3εt

The particle relaxation time τp is equal to τp )

(2(Fp /Ff) + 1)(2ξ)2 36νt

where Fp, Ff, and νf are the particle density, the fluid density, and the fluid kinematic viscosity, respectively. Finally, the parameter b is expressed as

5226

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

b)

3 2(Fp /Ff) + 1

Baldyga and Orciuch,13 while reviewing several models to compute the particle turbulent diffusivities, discussed the limitations of Tchen’s theory. Tchen’s model, although limited for practical applications, is used here for the sake of clarity, the focus being on demonstrating how the closure problem can be solved in the FCMOM approach. For the same reason, eq 13 is assumed to be valid for all times during the simulations, even though, rigorously, it holds only for short diffusion times. However, the methodology followed is general and can be applied to other (more realistic) diffusivity models. The values of the parameters in the simulations performed are as follows: the particle density is equal to 2000 kg/m3, while the air density and viscosity under atmospheric conditions are used for the fluid phase; TL is approximately equal to 0.1 s. The initial and boundary conditions are the same as in subsection 3.1, with ξmin(0,x) and ξmax(0,x) (initial minimum and maximum particle radius) equal to 10 and 50 µm, respectively, in the whole spatial domain. Figure 4 shows the simulation results at x/L ) 0.5 when the time is equal to (1/10)(L2/Dft). Again, the comparison between the numerical solution (obtained using the FCMOM) and the exact solution is excellent. Figure 4 shows also the solution for Dpt ) Dft (size-independent diffusivity), to emphasize the effect of the particle size on the size-dependent diffusion process. The two solutions (for the size-dependent case and for the sizeindependent case) will eventually (at longer times) approach each other. 4. Application to Particle Convection In this section, the FCMOM technique is applied to particle convection. The application considered is the gas-solid flow in a pipe. The initial and boundary conditions are chosen such that the flow is axisymmetric and the problem is solved in a 2-D geometry (radial and cylindrical coordinates r, z). The FCMOM eqs 10 are coupled with the equations governing the gas-phase fluid dynamics (mass and momentum balances). The gas phase is in a laminar regime, and the particles move only due to the convection (particle diffusion is neglected).

Figure 4. PSD vs particle radius: comparison between the numerical solution by the FCMOM and the exact solution for pure particle diffusion (sizedependent diffusivity). The corresponding solution for Dpt ) Dft (sizeindependent diffusivity) is also plotted. The location is at x/L ) 0.5, and the time is (1/10)(L2/Dft). The numerical solution is obtained with M ) 11.

It is further assumed that the particles are light (low Stokes number) and in dilute conditions. Accordingly, the particle-phase velocity is assumed to be equal to the gas velocity. As already mentioned, this is not a restriction of the methodology presented in this paper. Such assumptionsthat the particle-phase velocity is equal to the gas velocityswas made because, in this way, it is straightforward to check whether the results predicted by the FCMOM (the particle size distribution) are correct or not. However, the methodology proposed is valid also when the particle-phase velocity is not equal to the gas-phase velocity; if so, the particle-phase (convective) velocity can be computed from another equation, i.e., the particle-phase momentum balance equation (Gidaspow15). In subsection 4.1 the gas-phase governing equations are introduced. In subsection 4.2 the FCMOM eqs 10 are simplified for the case of pure convection. In the two final subsections the simulation results are presented, first in conditions such that the gas phase is at steady state (subsection 4.3) and then for unsteady-state conditions (subsection 4.4). 4.1. Gas-Phase Governing Equations. The equations governing the gas-phase fluid dynamics (in isothermal conditions) are the mass and momentum balances: ∂Fgε + ∇ · (Fgεvg) ) 0 ∂t

(14)

∂Fgεvg + ∇ · (Fgεvgvg) ) -∇p - ∇ · τ + Fgεg + Fpg ∂t

(15)

Figure 3. PSD vs particle radius: comparison between the numerical solution by the FCMOM and the exact solution for pure particle diffusion (sizeindependent diffusivity). The location is at x/L ) 0.5 and the final time is L2/Dpt. The numerical solution is obtained with M ) 11.

In eqs 14 and 15, Fg, ε, vg, p, and τ are the gas density, the gas volume fraction, gas velocity, the gas pressure, and the stress tensor, respectively. The constant g is the gravitation acceleration, and Fpg represents the force (per unit volume) exerted by the particle phase on the gas phase (drag force and buoyancy). The gas phase is assumed incompressible and Newtonian; therefore, the stress tensor can be expressed using the following constitutive equation (Bird et al.16): 2 τ ) 2µgD + kg - µg tr(D)I 3

(

)

(16)

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

5227

where µg and kg are the shear and dilatational viscosities of the gas phase, I is the identity tensor, and D is the gas-phase rate of deformation tensor: 1 D ) (∇vg + ∇vTg ) 2 In the simulations performed, the particle-phase flow is dilute and the effect of Fpg is assumed to be negligible. In the simulations, the properties of air were used for the gas-phase viscosities and density. 4.2. FCMOM Equations. The evolution equations of the dimensionless moments (10), for a particle convective flow, can be rewritten as ∂µi + ∇ · (µivp) ) -(MB + MBconv) ∂t

(17)

where, as mentioned, vp ) vg. The gas-phase volume fraction is related to the moments through the following equation: ε ) 1 - fsc

(ξmax - ξmin) µ0 2

Figure 5. PSD vs particle radius: evolution of the PSD at the center of the pipe outlet for particle convection. The numerical solution by the FCMOM is compared with the inlet PSD at the final time. The numerical solution is obtained with M ) 9.

(18)

4.3. Particle Size Distribution with Gas-Phase SteadyState Conditions. In all the simulations performed for pure convection, both the gas-phase governing equations and the FCMOM equations were solved. However, in this subsection, the initial and boundary conditions are such that the flow field of the gas phase is constant and homogeneous. The steady-state condition is not imposed, but it results directly from the simulations. A short pipe with radius R ) 0.05 m and length L ) 0.3 m is used in the simulations. The initial condition of the gas-phase axial and radial velocity components is

[ ( Rr ) ]

Vg,z(0, r, z) ) Vg,max 1 -

2

Vg,r(0, r, z) ) 0

(19)

i.e., a uniform laminar fully developed flow profile is prescribed in each location. The same profile is prescribed at the inlet

[ ( Rr ) ]

Vg,z(t, r, 0) ) Vg,max 1 -

2

Vg,r(t, r, 0) ) 0

(20)

At the outlet, Vg,r(t,r,L) ) 0 and a Neumann boundary condition is prescribed for the axial component of the gas phase; i.e., its derivative is equal to zero. The simulations were performed using Vg,max ) 0.2 m/s. Two sets of simulations were performed. In the first set, the PSD is initially uniform and equal to f(ξ, 0, r, z) ) a

{

}

q [ξ - ξmin(0, r, z)] × [ξmax(0, r, z) - ξmin(0, r, z)] [ξmax(0, r, z) - ξ] [ξmax(0, r, z) - ξmin(0, r, z)]

{

Figure 6. PSD vs particle radius: PSD in different locations of the pipe outlet for particle convection at the final time. The numerical solution by the FCMOM at the final time close to the pipe wall is compared with the inlet PSD at the initial time. The numerical solution is obtained with M ) 9.

}

s

(21)

with q ) 2 and s ) 4; additionally, ξmin(0,r,z) ) 10 µm and ξmax(0,r,z) ) 25 µm.

At the inlet, a Dirichlet boundary condition is prescribed; i.e., the PSD is still given by eq 12, with s ) 4. However, the parameter q varies linearly from the value of q ) 2 at the initial time to q ) 8 at t1 ) 5 s and ξmax(t,r,0) varies linearly from the value of 25 µm at the initial time to 50 µm at t1. Afterward the inlet condition of the PSD remains constant. The results of the first set of simulations are shown in Figures 5 and 6, obtained with 9 moments (M ) 9). In Figure 5, the PSD at the outlet (z ) L) in the pipe center (r ) 0) is shown at the final time tfin ) t1 + L/Vg,max. The final PSD is correctly predicted by the FCMOM. In fact, the PSD which is at the inlet at time t1, at the final time tfin reaches, by convection, the pipe outlet. It is also noteworthy that the PSD at the outlet is different for different values of the radial coordinate r. In fact, because of the axial velocity radial profile, the inlet PSD is convected with different velocities for different r. As a result, the final outlet PSD near the wall has changed little from the initial inlet

5228

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

Figure 8. Spatial (r-z) distribution of the gas velocity radial component Vg,r/Vg,max in the pipe after 3.5 s. Figure 7. PSD vs particle radius: stationary PSD radial profile for particle convection in steady-state conditions (for the gas phase). The inlet PSD is equal to the initial condition. The numerical solution of the PSD inside the pipe (at r ) 0.04 m, z ) 0.05 m or at r ) 0 m, z ) 0.05 m) by the FCMOM is compared with the inlet PSD (r ) 0.04 m, z ) 0 m or at r ) 0 m, z ) 0 m) at the final time. The numerical solution is obtained with M ) 9.

condition, while on the centerline it has fully responded to changes at the inlet and is equal to the final inlet PSD (Figure 6). Again, this is predicted correctly by the FCMOM numerical results. A second set of simulations was performed to investigate the PSD radial variations in steady (for the gas phase) and unsteady conditions (see next subsection 4.4). The presence of such variations is connected, in a pure convective flow, to the velocity radial component. The initial PSD is radially inhomogeneous and equal to f(ξ, 0, r, z) ) a

{

}

q(r) [ξ - ξmin(0, r, z)] × [ξmax(0, r, z) - ξmin(0, r, z)] [ξmax(0, r, z) - ξ] [ξmax(0, r, z) - ξmin(0, r, z)]

{

}

s(r)

(22)

with q(r) ) 2 + 6r/R, s(r) ) 4, ξmin(0,r,z) ) 10 µm, and ξmax(0,r,z) ) 25 µm. The PSD inlet condition is equal to the initial condition (22). The simulation results are shown in Figure 7. It is expected that, being in steady-state conditions (for the gas phase) with negligible radial components of the gas velocity, the PSD inside the pipe remains equal to the initial PSD, i.e., to the inlet PSD. This behavior is correctly predicted by the numerical results (Figure 7). Such behavior changes when radial components of the gas velocity are not negligible (subsection 4.4). 4.4. Particle Size Distribution in Unsteady Conditions. In this subsection, a different set of initial and boundary conditions are used (for the gas phase), and the gas phase is also in unsteady conditions. The gas is assumed initially at rest; i.e., both the radial and axial gas-phase velocity components are zero in the whole domain at the initial time. Then, an axial flow is prescribed at the inlet, namely, a parabolic laminar profile is imposed at the inlet, while the radial components at the inlet are still equal to zero. The inlet value of Vg,max switches from 0 m/s (at the initial time) to 0.2 m/s. The boundary conditions at the outlet are as in previous subsection 4.3; i.e., Vg,r(t,r,L) ) 0 and ∂Vg,z(t,r,L)/∂z ) 0. Therefore, the gas phase is initially at rest and eventually will reach a final stationary-state condition defined by the final

Figure 9. PSD vs particle radius: PSD in the pipe (r ) 0.032 m, z ) 0.075 m) for particle convection in unsteady-state conditions after 2.75 s. The PSD is not equal to the inlet PSD with the same r; instead it is similar to the inlet PSD at r ) 0.024 m, i.e., closer to the pipe center. The numerical solution is obtained with M ) 9.

inlet laminar profile (with Vg,max ) 0.2 m/s), in which again the velocity radial components are zero. However, during the transient, the flow is developing and the radial components of the gas velocity are not negligible. Figure 8 shows the spatial (r-z) distribution of the gas velocity radial component Vg,r/Vg,max after 3.5 s. The PSD initial and inlet conditions are the same as in previous subsection 4.3. However, because of the gas convective radial movement during the transient, the PSD inside the pipe in a specific point r0-z0 is not determined, as in the steadystate simulations, only by the inlet condition at r0-z ) 0. In fact, the PSD inside the pipe in a specific point r0-z0 is the PSD convected by an element of fluid whose position at the inlet may be different than r0; specifically, in the simulations performed, the gas velocity radial components are positive (particles move from the pipe center toward the pipe wall), and therefore, the PSD inside the pipe in a specific point r0-z0 is the PSD convected by an element of fluid whose starting position at the inlet was closer to the pipe center. This behavior is predicted by the FCMOM numerical results and is shown in Figure 9.

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

5. Conclusions The FCMOM is a numerical technique to solve monovariate and bivariate population balance equations (Strumendo and Arastoopour7,8). The FCMOM is computationally efficient, provides accurate values both for the moments and for the reconstructed particle size distribution (or bivariate particle distribution function), and, in the FCMOM, the domains of the internal variables are clearly defined within boundaries, thus assuring that the values assumed by the internal variables are always realistic. The FCMOM was previously formulated and validated for systems in homogeneous conditions (Strumendo and Arastoopour7,8). In this work the FCMOM was extended to monovariate inhomogeneous systems, thus providing the theoretical basis to use the FCMOM in CFD applications. In extending the FCMOM to inhomogeneous systems, the assumption that the particle-phase convective velocity is independent of the internal variables (particle size) was retained. Therefore, the formulation proposed is valid when the particle-phase convective velocity is assumed to be a function only of the external variables (time and space). The FCMOM methodology was applied to a general monovariate PBE, neglecting the source terms (due to particle growth, nucleation, aggregation, and breakage, etc.). The treatment of the source terms in the FCMOM framework was already considered by Strumendo and Arastoopour.7 The moments evolution equations were derived from the general formulation of the PBE; for inhomogeneous systems, the moving boundary conditions of the minimum and maximum particle size can be obtained by evaluating the spatial derivatives of the minimum and maximum particle size and the particle velocities. The method was illustrated by applications to particle diffusion and particle convection. In the case of particle convection, a gas-solid dilute flow in a pipe was simulated and the FCMOM equations were coupled with the governing equations (mass and momentum balances) of the gas phase. In all of the applications considered, the numerical solutions obtained by the FCMOM either converged to the exact solution or (when an exact solution was not available) provided results in agreement with the expectations. Nomenclature a ) constant in the initial particle size distribution function B ) birth of particles (per unit time, volume, particle size) c ) coefficients in the size distribution function series expansion (Legendre case) D ) death of particles (per unit time, volume, particle size) Dpt ) particle turbulent diffusivity Dft ) fluid turbulent diffusivity Dg ) gas-phase rate of deformation tensor Fpg ) force (per unit volume) exerted by the particle phase on the gas phase f ) particle size distribution function G ) particle growth/dissolution rate I ) identity tensor kg ) gas-phase dilatational viscosity kt ) turbulent kinetic energy L ) interval in a 1-D Cartesian domain or pipe length M ) number of terms in the truncated series expansion approximating the particle size distribution and number of moments in a simulation m ) summation index n ) summation index N ) number of particles per unit volume

5229

P ) Legendre polynomials p ) gas pressure q ) parameter in the initial particle size distribution function R ) pipe radius r, z ) cylindrical coordinates (2-D) rp ) particle (or droplet) radius s ) parameter in the initial particle size distribution function t ) time TL ) Lagrangian integral time scale for the fluid motion u ) velocity V ) velocity (convective) x ) vector of Cartesian coordinates x, y, z ) Cartesian coordinates Greek Letters ε ) gas volume fraction εt ) turbulent kinetic energy dissipation rate µ ) moment of the particle size distribution function µg ) gas-phase shear viscosity ν ) summation index νf ) fluid kinematic viscosity F ) density τp ) particle relaxation time τ ) stress tensor φ ) trial functions (Legendre case) Subscripts and Superscripts f ) fluid fin ) final value g ) gas i ) moment order in ) initial value max ) maximum value min ) minimum value p ) particle sc ) scale factor t ) turbulent -1, +1 ) dimensionless radius equal to -1 and to +1, respectively ) dimensionless quantities ′ ) notation for the function of the dimensionless internal variable

Literature Cited (1) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Systems in Engineering; Academic Press: London, 2000. (2) McGraw, R. Description of Aerosol Dynamics by the Quadrature Method of Moments. Aerosol Sci. Technol. 1997, 27, 255. (3) Barrett, J. C.; Webb, N. A. A Comparison of Some Approximate Methods for Solving the Aerosol General Dynamic Equation. J. Aerosol Sci. 1998, 29, 31. (4) Marchisio, D. L.; Pikturna, J. T.; Fox, R. O.; Vigil, R. D. Quadrature Method of Moments for Population-Balance Equations. AIChE J. 2003, 49, 1266. (5) Marchisio, D. L.; Vigil, R. D.; Fox, R. O. Quadrature Method of Moments for Aggregation-Breakage Processes. J. Colloid Interface Sci. 2003, 258, 322. (6) Marchisio, D. L.; Fox, R. O. Solution of Population Balance Equations Using the Direct Quadrature Method of Moments. J. Aerosol Sci. 2005, 36, 43. (7) Strumendo, M.; Arastoopour, H. Solution of PBE by MOM in Finite Size Domains. Chem. Eng. Sci. 2008, 63/10, 2624. (8) Strumendo, M.; Arastoopour, H. Solution of Bivariate Population Balance Equations Using the Finite Size Domain Complete Set of Trial Functions Method of Moments (FCMOM). Ind. Eng. Chem. Res. 2009, 48, 262. (9) Hulburt, H. M.; Katz, S. Some Problems in Particle Technology. Chem. Eng. Sci. 1964, 19, 555. (10) Fan, R.; Marchisio, D. L.; Fox, R. O. Application of the Direct Quadrature Method Of Moments to Polydisperse Gas-Solid Fluidized Beds. Powder Technol. 2004, 139, 7.

5230

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

(11) Fan, R.; Fox, R. O. Segregation in Polydisperse Fluidized Beds: Validation of a Multi-Fluid Model. Chem. Eng. Sci. 2008, 63, 272. (12) Courant, R.; Hilbert, D. Methods of Mathematical Physics; Interscience: New York, 1953. (13) Baldyga, J.; Orciuch, W. Some Hydrodynamic Aspects of Precipitation. Powder Technol. 2001, 121, 9. (14) Hinze, J. O. Turbulence; McGraw-Hill: New York, 1975. (15) Gidaspow, D. Multiphase Flow and Fluidization; Academic Press: San Diego, CA, 1994.

(16) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena; John Wiley and Sons: New York, 2002.

ReceiVed for reView September 9, 2009 ReVised manuscript receiVed February 15, 2010 Accepted February 22, 2010 IE901407X