Direct Quadrature Method of Moments for the Mixing of Inert

Ind. Eng. Chem. Res. 2010, 49, 5141–5152

5141

Direct Quadrature Method of Moments for the Mixing of Inert Polydisperse Fluidized Powders and the Role of Numerical Diffusion Luca Mazzei,*,† Daniele L. Marchisio,‡ and Paola Lettieri† Department of Chemical Engineering, UniVersity College London, Torrington Place, London WC1E 7JE, U.K., and Dipartimento di Scienza dei Materiali e Ingegneria Chimica, Politecnico di Torino, 10129, Torino, Italy

Computational fluid dynamics (CFD) is extensively employed to investigate dense fluidized suspensions. Most mathematical models assume that the powder is monodisperse or is formed by few solid phases of particles with constant size. Real powders, nevertheless, are polydisperse, with their particle size distribution continuously changing in time and space. To account for this important feature, models have to include a population balance equation (PBE), which needs to be solved along with the customary fluid dynamic transport equations. The recently developed direct quadrature method of moments (DQMOM) permits solving PBEs in commercial CFD codes at relatively low computational cost. This technique, nevertheless, still needs testing in the context of dense multiphase flows. In this work we implement DQMOM within the CFD code Fluent to study the mixing of two polydisperse fluidized suspensions initially segregated. Each node of the quadrature represents a distinct solid phase advected with its own velocity. Simulating this apparently simple system highlights a problem related to the numerical solution of the DQMOM transport equations: these do not feature diffusive terms, but the numerical diffusion generated by the finite-volume integration method alters the model predictions, leading to wrong results. To solve this, the PBE needs to account for diffusion: this yields source terms in the transport equations of the quadrature weights and nodes that ensure the latter are correctly predicted. 1. Introduction Fluidization is a well-established technology used in many industrial processes, such as coal combustion, biomass gasification, waste disposal, sulfide ores roasting, and food processing. Fluidized bed reactors are attractive for treating fluid-solid systems because they maximize the contact area between the phases and guarantee excellent heat and mass transfer rates. To design them, engineers rely on experimental correlations and pilot plants. Nevertheless, since these correlations are valid only for the specific units investigated, they cannot help us improve design and performance; for instance, they can tell us nothing about the effect of changing size or geometry, introducing internals, such as baffles or heat exchanger tubes, or repartitioning the feeds over more entry points. Pilot plants, on the other hand, are costly and time-consuming; furthermore, they do not always lead to adequate scale-up. For these reasons, thanks also to the high-speed computers and advanced computational fluid dynamics (CFD) codes available nowadays, the modeling and numerical simulation of fluidized mixtures have rapidly gained importance, since they could reduce the need for expensive pilot plants, improve process performance, and shorten time to market. Even the most advanced models, however, usually do not take into account polydispersity, neglecting that the disperse phase is formed by elements with varying properties, such as particle size. They instead assume that the powder is constituted of classes of particles with equal and constant diameters. Often the system is supposed to be monodisperse,1-10 even if more recently some research groups have performed computational studies of bidisperse mixtures, still characterized by particles of constant diameter, trying to predict particle mixing and segregation.11-19 * To whom correspondence should be addressed. Tel.: +44 (0)20 7679 7868. Fax: +44 (0)20 7383 2348. E-mail: [email protected]. † University College London. ‡ Politecnico di Torino.

In some cases even more than two particle phases have been considered, as in the work of Mathiesen et al.,20,21 where to describe more realistically the particle size distribution in circulating fluidized beds the authors accounted for three solid phases. Real multiphase systems are instead characterized by wide particle size distributions (PSDs) that change continuously owing to fluid-particle and particle-particle interactions. Predicting their evolution is paramount if we want to describe realistically fluidized suspensions, since the PSDs affect product quality, mixing and segregation patterns, chemical reaction rates, heat transfer rates, loss of unreacted material, and many other parameters. The constant particle size assumption limits the model flexibility: solids can mix and segregate, but variations in their diameters are not allowed for. In reality, conversely, particles can shrink, aggregate, break, and nucleate; accordingly, their size distribution changes continuously in time and space. Predicting this evolution, which depends on the local conditions wherein the system operates, is essential for a reliable description of the mixture behavior. We can achieve this by solving, along with, or in place of, the customary multifluid equations of conservation for mass and linear momentum,22-24 a population balance equation (PBE). Since the dimensionality of the latter differs in general from that of classical fluid dynamic equations, succeeding in this task is difficult. Few attempts to implement the PBE into multiphase CFD codes can be found in the literature. Olmos et al.25 used this method to simulate bubble columns. They considered 10 different size groups to represent the bubbles, but solved only the dynamic equation for the mixture to save computational time. All the discrete phases were therefore convected with the same velocity. Other researchers have applied similar strategies to investigate gas-liquid systems;26-28 dense gas-solid systems, conversely, have never been investigated, with the exception of the pioneering works of Fan et al.29 and Fan and Fox.30

10.1021/ie901116y  2010 American Chemical Society Published on Web 12/09/2009

5142

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

The existence and uniqueness of the solution of PBEs is a standard subject in mathematical textbooks;31 thus, we can be sure that PBEs can be, at least in principle, solved. There are, however, very few analytical solution strategies, and those that do exist are almost invariably for such simplified and idealized problems that cannot be adopted for real systems. Various techniques can solve PBEs numerically; for a comprehensive review we refer to Ramkrishna.32 Here we mention the class method, also known as the discretized population balance approach. This first integrates the PBE over subintervals in a partition of the internal particle state space, and then reduces the dimensionality of the PBE, transforming it into a set of ordinary differential equations expressing macroscopic balances for the number of particles within each interval. This set includes as many equations as the number of intervals, or granular classes, making up the partition. The major drawback of the method is that the overall mass and number of particles is conserved only in the limit of an infinite number of classes; also, the PBE solution strongly depends on the grid, and since the number of discretized equations rises with the number of classes, high computational times are necessary for integrating the PBE on a sufficiently fine grid. Another technique for solving PBEs is running Monte Carlo simulations. These are based on artificial realizations of the system behavior. Initially, a population of particles is in a state known in a statistical sense; while it evolves, the population undergoes deterministic (continuous) changes in particle state, described by ordinary differential equations, and random (discontinuous) changes in particle state, due to processes such as agglomeration, breakage, and nucleation, with specified probabilities. To create a realization of the system behavior, we artificially generate random variables that satisfy the specified probability laws of change. By generating numerous realizations, we can determine the expected behavior of the system by averaging all the sample paths. This technique is powerful, but it requires a lot of computing power; with the computational resources available today it is still not practical for real applications, as demonstrated by Zucca et al.,33 who ran Monte Carlo simulations on a limited number of compartments to mimic a spatially inhomogeneous system. Often engineers are only interested in some integral properties of the density function which describes the particle population. Such properties, called moments, may be important because they control the product quality or because they are easy to measure and monitor. The idea behind the method of moments is to derive transport equations for the moments of interest by integrating out all the internal coordinates from the PBE.34 The method is attractive, for the number of equations to be solved is small; however, the moment transport equations are unclosed, since for any given set of moments that the modeler wishes to track, the equations usually involve also higher-order moments external to the set.35 Because of this, the method has been scarcely applied. The quadrature method of moments (QMOM) and the direct quadrature method of moments (DQMOM), which approximate the system density function using quadrature formulas, overcome this problem. Turning integrals of the density function into summations, the quadrature formula eliminates the problem of closure.36 This is the common element that QMOM and DQMOM share. The difference is in how the methods compute the nodes and weights of the quadrature formula. Forcing these to agree with a set of independent lower-order moments,37 QMOM tracks the moments by integrating their transport equations and back-calculates nodes and weights. DQMOM,

conversely, directly tracks the latter, solving transport equations that govern their evolution. For monovariate distributions, i.e., distributions with only one internal coordinate, to back-calculate the quadrature nodes and weights from the moments of the density function we can adopt the product-difference algorithm of Gordon,38 which simply requires finding the eigenvalues of a real symmetric tridiagonal matrix. However, being based on the theory of canonical moments,39 the algorithm cannot be applied when a higher number of internal coordinates are present. The quadrature approximation must then be determined by resolving a nonlinear algebraic system at each grid point of the computational domain and at each time step. This renders QMOM unattractive, for the method loses much of its simplicity and efficiency. Tracking directly the nodes and weights, DQMOM is more immediate and computationally cheaper. Another limitation of QMOM is that the model uses a phase-average velocity of the solid phases (which coincide with the quadrature nodes) to solve the moment transport equations. Thus, all phases share the same velocity. For multiphase flows this is a serious limit that should be overcome, if we want to render QMOM able to handle cases where each phase is advected with its own velocity. In this work, to simulate the mixing of two inert polydisperse powders initially segregated, we implement a new version of DQMOM, based on a volume rather than on a number density function, into the multifluid model of the CFD code Fluent. Choosing this density function is convenient because the transport equations deal directly with volume fractions instead of number densities. The particles neither react nor agglomerate nor break; therefore, the particle size distribution changes only because the powders mix. This is a relatively simple problem, but its very simplicity is key to testing the method, understanding it better, and highlighting possible issues or limitations. We believe that before tackling more complicated problems, involving continuous and discontinuous changes in the particle internal coordinates, this analysis is crucial. The article is organized as follows. We first introduce the problem that we are going to investigate. Next, we describe experimental materials, apparatus, methodology, and findings. We then present the mathematical model and the numerics. Finally, we discuss the predictions of the CFD simulations; as we shall see, these are inaccurate, with the model failing to correctly predict the evolution of the quadrature nodes. This will make us critically reexamine the DQMOM equations, letting us appreciate their limitations and the critical role of numerical diffusion. We then revise the model, accounting for diffusion, and show how the new numerical results compare with the experimental findings. 2. Description of the Problem Consider a fixed bed made up of two superposed layers of polydisperse particles of equal density. Referred to as powders A and B, respectively, lower and upper layers differ only in particle size distribution, with the one on top having greater mean particle size. If the PSDs are reasonably narrow, the materials have well-defined minimum fluidization velocities; we denote the latter by uA and uB, where uB > uA. The dynamics of the fluidized suspension depends on the ratio u/uB, where u is the superficial fluid velocity. If u/uB > 1, the two powders thoroughly mix; at lower ratios, powder B partly mixes with powder A and partly segregates toward the bed bottom. Figure 1 portrays these two conditions. At high fluid fluxes the system reaches a configuration where only one powder, perfectly homogeneous, is present. At low fluid fluxes, conversely, two


5143

Figure 1. Different fluidization dynamics depending on the superficial fluid flux.

layers coexist: the lower contains the segregated biggest particles of powder B, while the upper is a mix of powder A and the remaining particles of powder B. Here we investigate the simplest case: perfect mixing. Choosing a superficial fluid velocity greater than uB, we thoroughly mix powders A and B. The resulting powder, referred to in Figure 1 as powder C, has a PSD that combines the two original ones. Adopting the direct quadrature method of moments, we propose to predict the new distribution and check that it reasonably agrees with the experimental one. 3. Experimental Materials, Apparatus, Methodology, and Results The powders are ballotini with density of 2500 kg/m3. Figure 2 reports their density and cumulative PSDs obtained by sieving. Using these data, we computed the Sauter mean diameter of the distributions and their coefficients of variation.33,40 For powders A and B these are equal to 88 µm, 273 µm, 0.16, and 0.20, respectively. According to the classification of Geldart,40 both powders have fairly narrow distributions. The experimental minimum fluidization velocities uA and uB are roughly equal to 1.00 and 6.40 cm/s, respectively. The experimental apparatus consisted of a nearly twodimensional (2D) plexiglass rectangular column, 600 mm high, 350 mm wide, and 10 mm deep, with a 3.5 mm thick permeable sintered bronze rectangular plate as distributor. Nitrogen was fed via two flow meters. Pressure taps were installed 100 mm apart along the height of the bed; from these we collected pressure readings via a digitron electronic manometer. A system of two interlocked on/off valves operated simultaneously was installed on the rig to evacuate instantaneously the fluidizing gas during the bed freeze tests (see below) that we ran to analyze the mixing and segregation in the bed. Figure 3 shows a photograph of the setup and a schematic representation of the rig. As shown in Figure 1, in the experiment we filled half the column with powder, forming two layers of equal height (15 cm): the lower of powder A and the upper of powder B. We placed the coarser powder on top to enhance the fluidization dynamics and verify if segregation occurs. At the high fluid fluxes used in this study the two powders should mix almost perfectly, turning into a uniform suspension; however, if any segregation should take place, this is because the coarser powder naturally tends to sink toward the bottom of the bed. Hence, if we had placed powder B below powder A, we would have

Figure 2. Experimental normal and cumulative particle size distributions for powders A and B.

hindered this motion. After loading the powders, we fed nitrogen for 10 min, a time sufficient to reach pseudostationary conditions (as we checked by observing the bed). To guarantee vigorous mixing, we set the superficial velocity of the fluid to 15 cm/s, which results in a ratio u/uB equal to 2.35. Successively, we performed the so-called bed freeze test: by operating the interlocked on/off valves shown in Figure 3, we abruptly cut off the gas supply to the bed and vented the gas in the windbox

5144


Figure 4. Experimental particle size distributions of the lowest (near the distributor plate) and highest (near the bed freeboard) layers of the bed, after fluidization and collapse. Figure 3. Photograph and schematic representation of the experimental apparatus. 1, Nitrogen tanks; 2, oil filter; 3, flow meters; 4, pressure taps; 5, fluidized bed; 6, electronic manometer; 7, on/off valve control switch. X, Freeboard; Y, fluidized bed; Z, windbox.

of the vessel to the atmosphere. With this configuration most of the interstitial fluid in the powder initially escapes through the vent, finding there a lower fluid dynamic resistance. Then, as the particles accumulate on the distributor plate and the pressure drop across the settled bed becomes greater than that in the upper bed, the gas escapes from the top of the bed. This process is very quick, so the powder retains its instantaneous distribution, as if it were frozen. We then split the resting bed into five layers of equal height (6 cm), collected each layer by means of a sampling probe, and finally sieved them to obtain their PSDs. We ran the experiment three times. As we expected, the two powders mixed almost perfectly, with the particle size distributions in the five layers being nearly identical. Figure 4 reports the PSDs of the most significant layers: the lowest, which directly lies on the distributor, and the highest, which separates the bed from the freeboard. The new PSDs are identical and seem to be obtained by juxtaposing the two original distributions reported in Figure 2A,C. This indicates very good mixing. 4. Multiphase Fluid Dynamic Model The particle population that we consider is characterized by diameter and velocity; there are two internal coordinates, one scalar and one vectorial, and the internal state space is fourdimensional. To describe the population of particles, we introduce a volume density function (VDF); denoted by fV, this is defined so that fV(ξ, W, x, t) dξ dW dx represents the expected volume of particles contained at time t in the physical volume dx around x with size ξ in the range dξ and velocity W in the

range dW. We shall denote the domains of variation of ξ and W by Ωξ and ΩV, respectively. To determine fV, we should solve a population balance equation written in the four-dimensional internal state space mentioned above. Since neither heterogeneous reactions nor particle attrition occurs, ξ does not vary continuously and the particles have zero velocity in size space; thus, the PBE is ∂fV + ∇x · (fVW) + ∇V · (fVW˙ ) ) hV ∂t

(4.1)

where W˙ is the particle acceleration and hV is a source term that accounts for discontinuous jumps in particle state space. Instead of solving directly eq 4.1, we approximate the VDF using a quadrature formula, which expresses fV as a summation of ν Dirac delta functions, writing ν

fV(ξ, W, x, t) ≈

∑ φ (x, t) δ[ξ - ξ (x, t)] δ[W - W (x, t)] i

i

i

i)1

(4.2) where ξi(x, t) and Wi(x, t) are the ith quadrature nodes and φi(x, t) is the ith quadrature weight. Equation 4.2 states that the particle population is represented by ν classes, the ith of which has volume fraction φi(x, t), diameter ξi(x, t), and velocity Wi(x, t). Therefore, the problem reduces to predicting the evolution in time and space of these 3ν functions. To determine their equations of change, we should insert eq 4.2 into the PBE and then apply suitable moment transforms.41 The equations for the φi(x, t) would resemble continuity equations, while those for the Wi(x, t) would resemble dynamical equations. Here, however, for the velocities Wi(x, t) we adopt a different technique: we suppose that the customary averaged dynamical equations for multiphase flows hold (refer to section 4.2) and that the velocity


fields ui(x, t) that they provide coincide with the quadrature nodes Wi(x, t). Thus, we do not resort to moment transforms to obtain their equations of change, and only need to derive the transport equations for the quadrature weights φi(x, t) and nodes ξi(x, t). 4.1. DQMOM Transport Equations. To derive the DQMOM transport equations, we start by integrating out the coordinate W from eq 4.1; this results in the reduced population balance equation: ∂ jf V + ∇x · (fjV〈W|ξ〉) ) hjV ∂t

(4.3)

Equation 4.8 expresses the PBE for a monovariate population of particles whose velocities are conditioned on the scalar property ξ and whose VDF fulfills eq 4.6. The unknowns are the functions φi(x, t) and ξi(x, t), which denote the weights and nodes of the quadrature approximation respectively. Their evolution in time and physical space is governed by eqs 4.9; these define the functions ciφ(x, t) and ciφξ(x, t), but can be regarded as transport equations for the quadrature weights and weighted nodes. The source terms ciφ(x, t) and ciφξ(x, t) are unknown, but we can determine them by computing 2ν moment transforms of the PBE. Given a function φ(ξ, x, t), the integral transform

where, by definition, it is jf V ≡

∫

f ΩV V

∫

jf V〈W|ξ〉 ≡

dW;

fW ΩV V

dW;

hjV ≡

∫

ΩV

(4.4) Here 〈W|ξ〉(ξ, x, t) is the mean velocity conditioned on the particle size ξ. Since 〈W|ξ〉 is size-dependent, in eq 4.3 no diffusive flux in physical space is present. This is because particles with different sizes are advected with different velocities. Spatial diffusion would arise if we replaced 〈W|ξ〉 with the mean velocity of the whole particle population, this being 〈〈W|ξ〉〉(x, t) ≡

1 φ(x, t)

∫

j

f (ξ, x, t) Ωξ V

〈W|ξ〉(ξ, x, t) dξ (4.5)

where φ(x, t) is the overall particle volume fraction. For details we refer to Marchisio and Fox.41 If we use eq 4.2 to approximate the original density function fV, we find jf V(ξ, x, t) ≈

ν

∑ φ (x, t) δ[ξ - ξ (x, t)] i

i

and

i)1

〈W|ξ〉[ξi(x, t), x, t] ) Wi(x, t)(4.6) Hence, as expected, the particles belonging to the size class ξi move with velocity Wi(x, t), this coinciding with their conditional velocity 〈W|ξ〉(ξi, x, t). Note that because particles neither aggregate nor break, the size-dependent source term hjV(ξ, x, t) vanishes, and eq 4.3 becomes ∂ jf V + ∇x · (fjV〈W|ξ〉) ) 0 ∂t

i

φξ i

- ξicφi )δ′(ξ - ξi)] ) 0

(4.10)

ν

(1 - k)

∑c ξ

φ k i i

i)1

ν

+k

∑c

φξ k-1 i ξi

)0

(4.11)

i)1

This is an algebraic equation in the 2ν unknowns ciφ and ciφξ. To find these functions, we need 2ν independent equations. These can be obtained by writing eq 4.11 for 2ν different values of the parameter k; since the equations are in general coupled, a linear algebraic system of order 2ν must be solved. In this case, since the system is homogeneous, the solution is trivial: the source terms ciφ and ciφξ are zero, regardless of which set of moments we employ. This is because particles do not react, wear, nucleate, aggregate, or break. The selection of the moments that we wish to preserve comes in only when we initialize the functions φi(x, t) and ξi(x, t), for their initial values depend on the moments we choose (refer to section 6). 4.2. Multifluid Dynamical Equations. We assume that the dynamical equations for the fluid and particle phases (the latter really representing the quadrature nodes) are the customary multifluid equations of multiphase systems, obtained by mathematical averaging.42 Thus, for the fluid we have Fe

[ ∂t∂ (εu ) + ∇ · (εu u )] ) ∇ · S - ∑ n f + εF g e

e e

e

i i

e

i)1

(4.12) where Fe and ε are its density and volume fraction respectively, ue is its averaged velocity, and Se is its effective stress tensor. Moreover, ni is the number density of particle phase i and fi is the force exerted by the fluid on a single particle of the ith phase. Finally, g is the gravitational field. We do not need a transport equation for ε because ε ) 1 - φ, where φ is the sum of all the quadrature weights. The dynamical equation for the ith particle phase takes the form Fi

[

∂ (φ u ) + ∇ · (φiuiui) ) ∇ · Si + ni fi + ∂t i i

]

ν

∑n f

i ik

+

k)1

φiFi g(4.13)

(4.8)

i)1

where δ′(ξ - ξi) is the derivative of δ(ξ - ξi) with respect to ξ - ξi. Also, by definition, we have ∂φi + ∇x · (φiWi) ≡ cφi (x, t); ∂t

φ(ξ, x, t)ξk dξ

defines the moment of order k with respect to the internal coordinate ξ of the function φ(ξ, x, t). If we apply this transform to eq 4.8, we obtain

(4.7)

ν

φ i

Ωξ

ν

This new population balance equation governs the evolution of the monovariate VDF jfV(ξ, x, t). We can use it to find the 2ν functions φi(x, t) and ξi(x, t), but the equation can no longer provide any information about the velocities Wi(x, t). To compute these, we resort to the averaged dynamical equations of multiphase flows (refer to section 4.2), assuming as we said that Wi(x, t) ) ui(x, t). With this hybrid approach the particle velocity is no longer an internal coordinate and the VDF becomes monovariate, but additional equations are necessary in addition to the PBE, this being left only with the task of governing the evolution of the PSD. Now, inserting eq 4.6 into eq 4.7 yields

∑ [c δ(ξ - ξ ) - (c

∫

φ(ξ, x, t) f Mk(φ) (x, t) ≡ hV dW

5145

∂ (φ ξ ) + ∇x · (φiξiWi) ≡ ∂t i i cφξ i (x, t)(4.9)

where Fi is the solid density of phase i (which is the same for all the classes), Si is the effective stress tensor of phase i, fik is the force exerted by phase k on a single particle of phase i, and ui(x, t) is the velocity of phase i, which we assume to be equal to the quadrature node Wi(x, t). These dynamical equations are unclosed; to render them solvable, we need to relate the effective stress tensors Se and Si and the interphase forces fi and fik to the averaged velocity fields.

5146


4.2.1. Effective Stress. We close the effective stress tensors using customary Newtonian constitutive equations;43,44 accordingly, we write 2 Se ) -peI + 2µeDe + κe - µe tr DeI; 3

(

)

2 Si ) -piI + 2µiDi + κi - µi tr DiI (4.14) 3

(

)

where pe, pi, µe, µi, κe, and κi are the averaged pressures, shear viscosities, and dilatational viscosities of the fluid and particle phases, respectively; moreover, I is the identity tensor, while De and Di are the rate of deformation (or strain) tensors, which are defined as 1 1 De ≡ (∇ue + ∇ueT); Di ≡ (∇ui + ∇uiT) (4.15) 2 2 We regard the fluid as incompressible and do not specify its pressure constitutively; furthermore, we assume that µe is constant and neglect κe. The solid phases are characterized in the viscous flow regime by a viscous solid pressure piv, shear viscosity µvi , and dilatational viscosity κvi , and in the plastic flow regime by a plastic pressure pip and shear viscosity µip. In the viscous regime the generic property fi coincides with fvi , whereas in the plastic regime we assume fi ) fiv + fip. We express piv using the closure of Gidaspow:23

[

ν

pvi ) 1 + 2

∑ k)1

]

()

ξik 3 (1 + eik)φk gik φiFiΘi ξi

ξi + ξk (4.16) 2 Here Θi is the granular temperature of the ith phase, related to the kinetic energy of the fluctuating particle motion. eik is the restitution coefficient for the collisions between particles of phases i and k; for k ) i, eii coincides with the restitution coefficient ei of phase i. Finally, gik represents the radial distribution function that we obtain by combining the radial distribution functions gi and gk of the ith and kth particle phases, respectively. We close these functions as

[ ( )]

ξi ν φk φ + 12 k)1 ξk φmax

∑

1/3 -1

;

gik )

ξi gk + ξk gi ξi + ξ k

(4.17) φmax is the maximum solid packing. Still using Gidaspow’s closures, we have

( )

Θi 1/2 4 ; κvi ) φi2Fiξi gi(1 + ei) 3 π 10Fiξi√πΘi 2 4 3 µvi ) 1 + (1 + ei)φi gi + κvi (4.18) 96(1 + ei) gi 5 5

[

]

The granular temperatures are governed by balance equations for pseudointernal energies related to the particle velocity fluctuations.22,23 For fluidized mixtures these equations differ from the classical internal energy balance equation because of a sink term Sic representing losses of pseudointernal energy caused by inelastic collisions, a source term Gid representing the generation of particle velocity fluctuations by fluctuating fluid-particle forces, and a sink term Siv representing their dampening by the viscous resistance to particle motion. Thus, each balance equation is:

[ ∂t∂ (φ U ) + ∇ · (φ U u )] ) -∇ · q + S :∇u + G

Fi

i

i

i

i i

i

i

µpi )

pi sin ϑi 2√I2(Di)

pvi sin ϑi

)

2√I2(Di)

i

d i

- Svi - Sci

(4.19)

1 I2(Di) ≡ [(tr Di)2 - tr Di2] 2

;

(4.20) where ϑ is the angle of internal friction of the ith granular material and I2(Di) is the second invariant of Di. In eq 4.19 the higher viscosity generates higher dissipation of mechanical energy in pseudointernal energy, increasing the granular temperature and in turn the viscous solid pressure. 4.2.2. Fluid-Particle Interaction Forces. The fluid-particle interaction force consists of buoyancy and drag forces; thus, it is fi ) fis + fid. We define ni fis ≡ -φi∇pe and close ni fid as46 ni fid ≡ βi(ue - ui);

where ξik ≡

gi )

where Ui ≡ 3Θi/2 is the pseudointernal energy per particle unit mass and qi is the pseudothermal heat flux. The closure problem then requires finding constitutive expressions also for qi, Siv, Sic and Gdi . We do not report them here for briefness but refer to Gidaspow.23 To close the plastic solid stress, we use the so-called KTGFbased model. When φ exceeds the frictional threshold φf, we keep on using the viscous closure for the solid pressure, eq 4.16, but increase the solid shear viscosity by adding to the viscous contribution, eq 4.18, a frictional one whose expression was developed by Schaeffer:45

Fe |ue - ui |εφi -ψ(ε,Rei) 3 βi ) CD(Rei) ε 4 ξi (4.21)

where it is ψ(ε, Rei) ≡ -

ln φ(ε, Rei) ; ln ε

φ(ε, Rei) ≡

CD*(ε, Rei) 2(1-n) ε CD(Rei)

CD(Rei) ) (0.63 + 4.8Rei-1/2)2 ; CD*(Rei*) ) (0.63 + 4.8Rei*-1/2)2 Rei ≡

Fe ε|ue - ui |ξi ; µe

Rei*(ε, Rei) ≡

n(Rei*) )

Rei

; εn 4.8 + 2.4(0.175Rei*3/4)

(4.22) 1 + 0.175Rei*3/4 Here CD and CD* are drag coefficients evaluated with the correlation of Dallavalle,47 whereas n is the Richardson and Zaki48 coefficient evaluated with the correlation of Rowe.49 4.2.3. Particle-Particle Interaction Forces. We assume that the interaction force fik exchanged between particles of different phases includes only a draglike contribution, being thus proportional to the slip velocity between the phases. To close it, we adopt the constitutive equation of Syamlal:50 ni fik ≡ ζik(uk - ui); FiFkφiφk gik(ξi + ξk)2 3 π ζik ) (1 + eik) 1 + Fik |uk - ui | 4 4 Fξ3 + F ξ 3

(

)

i i

k k

(4.23) Here eik is set to 0.90, Fik is a coefficient of friction equal to 0.15, and gik is the radial distribution function defined in eq 4.17. 5. Numerical Schemes and Techniques To run the simulations, we employed the commercial CFD code Fluent 6.3. The governing and constitutive equations were implemented in the multifluid model of the package, which is


based on a Eulerian description of the flow. We treated the quadrature weighted nodes as user scalars and added their transport equations to the default equations of the code. We implemented the closure for the fluid-particle interaction force and postprocessed the simulation results by means of user functions (that is, routines and algorithms that run along with the CFD simulation). We used the pressure-based solver, which is recommended for low-speed incompressible flows. To convert scalar differential equations into algebraic equations which can be solved numerically, the code adopts a finite-volume discretization scheme. To ensure convergence, we discretized in space through a first-order upwind scheme, where cell-face quantities are determined by assuming that the cell-center values of any field variable represent cell averages that hold throughout the entire cells; thus, face quantities are identical to cell quantities and are set equal to the cell-center values in the upstream cells (relative to the direction of the normal velocity). Temporal discretization is first-order accurate and implicit. To couple pressure and velocity, we adopted the SIMPLE (Simultaneous Solution of Nonlinearly Coupled Equations) algorithm of Lo;51 the code does not allow any other coupling algorithms for Eulerian multiphase calculations. At each time step we employed a maximum of 150 iterations to compute the flow variables. By setting the tolerance to 10-5 we attained convergence within the iteration limit. The time step was set to 10-2 s, since we ascertained that shorter time steps yielded the same results. Under-relaxation factors of 0.20 were adopted for all the variables. 6. Boundary and Initial Conditions The computational grid (uniform, with square cells of 5 mm side) is two-dimensional; hence, front and back wall effects were neglected. On the left and right walls we used no-slip boundary conditions. At the bottom of the bed the inlet gas velocity was set to 15 cm/s. At the domain upper boundary the pressure was set to 105 Pa. On all boundaries the quadrature weights and weighted nodes fluxes were set to zero. To assign the initial conditions, we need to know the values that the 2ν quadrature weights and nodes φi(x, t0) and ξi(x, t0) have at the initial time t0. To this purpose, we must know 2ν independent moments of the VDF everywhere within the computational domain. Once we know them, we can use the quadrature formula of eq 4.4 to write for the generic kth moment: Mk(x, t0) ≡

∫

Ωξ

jf V(ξ, x, t0)ξk dξ≈

ν

∑ φ (x, t ) ξ (x, t ) i

0

k i

0

(6.1)

i)1

This is a nonlinear equation in the unknowns φi(x, t0) and ξi(x, t0). If we write eq 6.1 for 2ν values of k, we can then solve the resulting nonlinear algebraic system and consequently determine the values of φi(x, t0) and ξi(x, t0). As said, the solution depends on the moments that we use in the calculation: different choices initialize the problem differently and do not preserve the same properties of the VDF. In this work we employed a two-node quadrature approximation; as a consequence, we had two weights and two nodes to determine. Four unknowns require four moments. We decided to use the first four integer moments for two reasons. First, this choice turns eq 4.4 into a Gaussian quadrature, mathematically more accurate.41,52 Second, the selected moments relate to important properties of the distribution: the overall solid volume fraction, the mean particle size, and the VDF variance and skewness.34,40 Note that, since the quadrature approximation is Gaussian, we can determine weights and nodes very efficiently

5147

Table 1. Values of the VDF Moments and of the Quadrature Nodes and Weights Obtained from the Experimental PSDs Reported in Figure 2A,C Assuming a Void Fraction of 0.400 Moments of the Volume Density Function powder

M0

M1 [µm]

M2 [µm2]

M3 [µm3]

A B

0.600 0.600

5.45 × 101 1.70 × 102

5.06 × 103 4.98 × 104

4.82 × 105 1.52 × 107

Quadrature Nodes and Weights powder

ξ1 [µm]

φ1

ξ2 [µm]

φ2

A B

75 240

0.262 0.380

103 355

0.338 0.220

by adopting the product-difference algorithm of Gordon38 to solve the nonlinear algebraic system. In its initial state the bed is packed and constituted of two superposed layers; these are 15 mm high and together occupy half of the vessel. Since we know the experimental PSDs in the two layers (refer to Figure 2A,C), we can determine the moments Mk(x, t0); note that, the powders being well mixed, in each layer the moments do not depend on x. To calculate Mk, we use the relation n

Mk ≈ (1 - ε)

k+1 ω(di-1, di) dik+1 - di-1 k+1 i - di-1

∑d i)1

(6.2)

where n is the number of sieves used (10 in our case), di is the aperture of the ith sieve, and ω(di-1, di) is the mass fraction of powder in the size range (di-1, di). Equation 6.2 tells us that Mk is a function of ε; this is because, whereas the PSD refers to solid mass fractions on a void-free basis, the VDF accounts for voids and provides volume densities, that is, volumes of solid per unit volume of physical space. Table 1 reports the values of moments, nodes, and weights that we calculated by setting in eq 6.2 ε ) 0.40. With a twonode representation powders A and B become bidisperse: the first with particles of diameters 75 and 103 µm and volume fractions 0.262 and 0.338 respectively, whereas the second with particles of diameters 240 and 355 µm and volume fractions 0.380 and 0.220 respectively. We used these values to initialize the quadrature weights and weighted nodes. For visual clarity Figure 5 pictures the two-node representations of the PSDs; the arrows symbolize the Dirac delta functions of the quadrature formulas, their positions and heights indicating the quadrature nodes and weights on a void-free basis (that is, the mass fractions of the quadrature classes) respectively. For reference the figure also reports the experimental PSDs obtained by sieving. 7. Results and Discussion As previously said, in the experiments the powders perfectly mix. From the resulting PSD reported in Figure 4 we can calculate the first four integer moments of the VDF, using eq 6.2, and then back-calculate the quadrature weights and nodes, using the product-difference algorithm of Gordon.38 Table 2 shows the results, where we assumed a bed voidage equal to 0.40. The two-node representation again models the powder as a bidisperse suspension; the two sets of particles have diameters 105 and 304 µm and volume fractions 0.331 and 0.269 respectively. Note that, since the particles are inert and do not break, agglomerate, or nucleate, and since perfect mixing is attained, the nodes take on values which are intermediate between their original ones: from 75 and 240 µm the first node tends to 105 µm, whereas from 103 and 355 µm the second

5148


Figure 6. Two-node representation of the experimental PSD of powder C (obtained by mixing powders A and B). The arrows symbolize the Dirac delta functions of the quadrature formula, their positions and heights indicating respectively the quadrature nodes and weights on a void-free basis. For reference, the figure also reports the experimental PSD obtained by sieving.

Figure 5. Two-node representations of the experimental PSDs of powders A and B. The arrows symbolize the Dirac delta functions of the quadrature formulas, their positions and heights indicating respectively the quadrature nodes and weights on a void-free basis. For reference, the figure also reports the experimental PSDs obtained by sieving. Table 2. Values of the VDF Moments and of the Quadrature Nodes and Weights Obtained from the Experimental PSD Reported in Figure 4 Assuming a Void Fraction of 0.400 Moments of the Volume Density Function M0

M1 [µm]

M2 [µm2]

M3 [µm3]

0.600

1.16 × 102

2.85 × 104

7.95 × 106

Quadrature Nodes and Weights ξ1 [µm]

φ1

ξ2 [µm]

φ2

105

0.331

304

0.269

node tends to 304 µm. Figure 6 pictures the PSD two-node representation, also in this case superposed, for reference, to the experimental PSD obtained by sieving. The numerical simulation comprises two serial stages: a first where the suspension is fluidized by feeding gas at a superficial velocity of 15 cm/s, and a second where the gas supply is cut off and the system settles again into a packed bed. To let the simulation reach pseudostationary conditions, we simulated the fluidized bed for 10 s before cutting off the gas supply; the settling phase took instead about 2 s to complete. Figure 7 reports the profiles of the quadrature nodes and weights at the start of the simulation, at pseudo steady state, and at the end of the simulation. Phases 2 and 3 represent the first and second quadrature classes respectively, with the second corresponding to the greater particle size. At the high fluid flux used the bed dynamics is very fast; the volume fraction profiles are not uniform, for the system operates in the bubbling regime, but vigorous mixing occurs continuously. This was observed both in the experiment and in the simulation. After having let the bed settle, we divided it into

five layers and therein averaged the quadrature weights and nodes, mirroring what we did experimentally. Table 3 reports the values of the quadrature weights in the highest and lowest layers, showing also their averaged values over the entire bed. The lowest layer is more densely packed, having a voidage of 0.360, which is the maximum packing limit used in the simulation; the highest layer has instead a voidage of 0.440. To see if the powder had partly segregated, we had thus to refer to void-free mass fractions; these are also reported in Table 3 and are nearly identical in both layers. The top one is slightly richer in the smaller particles. This could be because the bed settles instantaneously in the experiments but takes 2 s in the simulation. During this time segregation might start. In any case the difference is negligible. Notice that in the simulation, after the gas supply is cut off, the interstitial gas leaves through the bed surface, since no vent is present, and this slows down the bed collapse. The quadrature nodes have uniform profiles, being equal to 165 and 190 µm everywhere in the bed. These values do not vary while the bed settles. Table 4 shows the experimental and numerical values of the quadrature nodes and weights: the weights are similar (error within (5%), but the nodes are quite different (error within (60%); qualitatively their values are acceptable, since ξ1 is between 75 and 240 µm and ξ2 is between 103 and 355 µm, but quantitatively they are unsatisfactory. Figure 8 helps to visualize the computational results; in the background the figure also reports the experimental PSD obtained by sieving. As we can see, the Dirac delta functions are not located as in Figure 6, even if the heights are roughly the same. This confirms that the weights are correctly predicted whereas the nodes are not. To explain this discrepancy, let us examine the DQMOM transport equations. In Eulerian form they are expressed by eqs 4.9; in Lagrangian form they become Diφi ) -φi∇ · Wi ; Dt

Diξi ) 0 where Dt

Di( · ) ∂( · ) ≡ + Dt ∂t Wi · ∇( · )(7.1)

The first equation tells us that in each differential material element the quadrature weights vary only if the divergences of their velocity fields are not zero. If the elements expand or contract, the quadrature weights change accordingly to preserve the phase volumes. The second equation tells us that along the stream lines the quadrature nodes remain constant to their initial


5149

Figure 7. Profiles of the quadrature weights and nodes at the start of the simulation, at pseudo steady state, and at the end of the simulation. Phases 2 and 3 represent the first and second quadrature classes, respectively. Table 3. Computational Quadrature Weights and Void-Free Mass Fractions in the Lowest and Highest Layers of the Bed and Averaged over the Entire Bed 2 s after the Feed Cutoff

highest layer lowest layer entire bed

ε

φ1

φ2

ω1

ω2

0.440 0.360 0.400

0.302 0.340 0.320

0.258 0.300 0.280

0.539 0.531 0.534

0.461 0.469 0.466

values. The motion rearranges the elements within the domain, causing what we usually call macromixing. Each element, however, evolves independently, without interacting with the other elements; accordingly, no micromixing takes place. This is what we would observe if we were able to solve exactly the equations reported above. This highlights a clear limitation of

Table 4. Experimental and Computational Values of the Quadrature Nodes and Weightsa

experimental computational

ε

φ1

φ2

ξ1 [µm]

ξ2 [µm]

0.400 0.400

0.331 0.320

0.269 0.280

105 165

304 190

a The numerical results derive from the original DQMOM transport equations, which do not account for diffusion and relevant source terms.

the model, which, treating the quadrature nodes as passive, nondiffusive scalars, cannot predict the microscopic mixing that really occurs. The scalar fields ξi(x, t) cannot become uniform and remain segregated, with the bulk motion only lowering the segregation length scale.

5150

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010 Table 5. Experimental and Computational Values of the Quadrature Nodes and Weightsa

experimental computational

ε

φ1

φ2

ξ1 [µm]

ξ2 [µm]

0.400 0.400

0.331 0.328

0.269 0.272

105 103

304 312

a The numerical results derive from the revised DQMOM transport equations, which account for diffusion and relevant source terms.

Here i and j denote the two quadrature classes, while the terms Ci and Cj are defined to be Ci ≡ Dxφi(∇xξi · ∇xξi);

Figure 8. Two-node representation of the computational PSD of powder C (obtained by mixing powders A and B) yielded by the original DQMOM transport equations, which do not account for diffusion and related source terms. The arrows symbolize the Dirac delta functions of the quadrature formula, their positions, and their heights indicating respectively the quadrature nodes and weights on a void-free basis. For reference, the figure also reports the experimental PSD obtained by sieving.

These considerations suggest that, at least in the transport equations of the quadrature nodes, a diffusive contribution should be present. Even so, we cannot motivate this on theoretical grounds. As already pointed out, each quadrature class moves with its own velocity, which is governed by its own dynamical equation. When this modeling approach is used, diffusion does not arise. If we considered a common velocity for all quadrature classes, diffusion would instead appear to make up for the difference between the actual velocity of each class and the mean velocity with which their advection is modeled. For more details about this point we refer to ref 44 and to the clear analysis of Curtiss and Bird.53 Notwithstanding, the numerical solution that we do obtain is quite different, as Figure 7 clearly shows: here micromixing does occur, since the quadrature nodes change tending toward a common value. This is caused by numerical diffusion: only diffusion consents micromixing, and since no diffusive terms appear in the transport equations, its origin must be numerical. The results are therefore incorrect because we are not solving the real equations of the model. The reader, however, might wonder why the numerical values of the micromixed nodes are wrong; numerical diffusion should simply mix the nodes, letting them evolve toward their correct mean values. To understand why this does not happen, reconsider the original PBE, eq 4.7; if we introduce diffusive transport, this becomes ∂ jf V + ∇x · (fjV〈W|ξ〉) - ∇x · (Dx∇xjf V) ) 0 ∂t

(7.2)

where Dx is a diffusion coefficient in physical space introduced by the numerics. With a technique similar to that presented in section 4.1, we can now derive new DQMOM transport equations; these are ∂φi + ∇x · (φiWi) - ∇x · (Dx∇xφi) ) cφi ∂t

(7.3)

∂ (φ ξ ) + ∇x · (φiξiWi) - ∇x · (Dx∇x(φiξi)) ) cφξ (7.4) i ∂t i i where the source terms on the right-hand side are equal to cφi )

6(Cj - Ci) (ξi - ξj)

2

;

cφξ i )

2Cj(ξj + 2ξi) - 2Ci(ξi + 2ξj) (ξi - ξj)2

(7.5)

Cj ≡ Dxφj(∇xξj · ∇xξj)

(7.6)

We then see that the diffusive flux ∇x · (Dx∇x fV) in the PBE not only generates diffusion in the transport equations of the quadrature weights and weighted nodes, but also generates source terms. In our physical problem, since neither continuous nor discontinuous processes change the particle diameters, ciφ and ciφξ should be zero, as they were in eq 4.9. However, the equations above reveal that in the presence of numerical diffusion they do not vanish, and relate to the diffusion coefficient Dx and the node gradients ∇xξi and ∇xξj. Thus, if we have diffusion, we also have generation. The source terms are essential, for they make up for the spurious terms, yielded by numerical diffusion, that appear in the quadrature transport equations. For further details, we refer to Fox54 and Marchisio and Fox.36 In Lagrangian form, the quadrature transport equations become Diφi ) -φi∇ · Wi + ∇x · (Dx∇xφi) + cφi Dt Diξi ) ∇x · (Dx∇xξi) + cξi ; Dt

cξi ≡

φ cφξ i - ξici φi

(7.7) (7.8)

The material elements are no longer segregated, and in addition a source term appears. In inhomogeneous systems neglecting the latter makes the nodes evolve toward incorrect values. When we solved eqs 4.9, even if diffusion was absent, the numerical scheme generated it; thus, albeit indirectly, it was taken into account. This, as pointed out, made it possible for the nodes to micromix. The source terms, though, were missing and this explains why the final values of the quadrature nodes were wrong. Equations 7.3 and 7.4 should yield the correct results, provided we estimate correctly the coefficient of numerical diffusion so that the source terms ciφ and ciφξ that we have to implement into the code are consistent with the diffusive fluxes ∇x · (Dx∇xφi) and ∇x · (Dx∇x(φiξi)) generated numerically. To estimate Dx, we used the relation reported by Ferziger and Peric:55 Dx )

uLc 2

(7.9)

where u is the velocity at which the property is convected and Lc is the length of the computational cell. This relation is valid only for the first-order upwind discretization scheme, which we adopted in the simulation. Taking as characteristic velocity 0.10 m/s, a value that has the same order of magnitude as the gas superficial velocity, and Lc equal to 5 mm, we obtain a diffusion coefficient of about 10-4 m2/s. If we solve the revised DQMOM transport equations using the value for Dx just computed, the fluidized suspension behaves qualitatively as shown in Figure 7; hence, we do not report again the snapshots of the computational profiles. The results this time agree with the experimental data; Table 5 reports them and, as we can see, the match is now more than


5151

could reasonably assume that both source terms and diffusive fluxes were based on the same coefficient. The results of the simulation are ξ1 ) 100 µm and ξ2 ) 321 µm, in good agreement with the experimental findings. 8. Conclusions

Figure 9. Two-node representation of the computational PSD of powder C (obtained by mixing powders A and B) yielded by the revised DQMOM transport equations, which account for diffusion and related source terms. The arrows symbolize the Dirac delta functions of the quadrature formula, their positions, and their heights indicating respectively the quadrature nodes and weights on a void-free basis. For reference, the figure also reports the experimental PSD obtained by sieving.

satisfactory. Note that the contour plots for the two nodes are still flat, showing that the source terms do not erase the effect of numerical diffusion, but only guarantee that the nodes are correctly predicted. Figure 9 visualizes the new computational results, showing in the background the experimental PSD. This time the Dirac delta functions are located as in Figure 6. If we observe the analytical expressions of the source terms ciφ and ciφξ, we see that they vanish when the gradients of the quadrature nodes are zero. This explains why in the recent work of Fan and Fox30 the problem that we detected passed unnoticed. Since they applied DQMOM to simulate how a perfectly mixed polydisperse fluidized powder segregates, not accounting for the source terms did not affect the solution because the node gradients and in turn ciφ and ciφξ were identically zero. In their simulations the quadrature nodes did not change from their initial values; segregation was reflected only by the uneven distribution of the solid volume fractions (quadrature weights) throughout the bed. Other cases where neglecting the source terms generated by numerical diffusion might not compromise the solution are when the particle sizes vary continuously and/or discontinuously because of heterogeneous reactions, attrition, nucleation, breakage, or aggregation. These and similar phenomena generate their own source terms ciφ and ciφξ in eqs 7.3 and 7.4, and if these terms outweigh those related to numerical diffusion, the resulting error might impact negligibly. To conclude the work, we conducted a sensitivity analysis on the coefficient of numerical diffusion. What happens if we do not estimate it correctly and fail to render the source terms on the right-hand side of eqs 7.3 and 7.4 consistent with the diffusive fluxes on the left-hand side? To answer this question, we ran two other simulations using in the source terms diffusion coefficients equal to 10-5 and 10-3 m2/s (note that the diffusive fluxes are generated by the code, so we did not have to implement them). The results that we found are ξ1 ) 78 µm and ξ2 ) 110 µm and ξ1 ) 47 µm and ξ2 ) 297 µm, respectively. As we expected, both simulations do not predict correctly the nodes. This confirms that diffusive fluxes and source terms must be consistent, that is, must both be based on the same diffusion coefficient. As a final check, we ran a simulation where we implemented both source terms and diffusive fluxes, using a diffusion coefficient equal to 10-2 m2/ s. Since this value is much bigger than the coefficient of numerical diffusion (which we estimated to be 10-4 m2/s), we

In this work, we presented a novel formulation of DQMOM based on a volume rather than on a number density function; this results in transport equations for the quadrature weights and nodes written in terms of volume fractions. We adopted the model to simulate the mixing of two inert polydisperse fluidized powders initially segregated. The equations did not predict correctly the evolution of the PSD. The problem has to do with the mathematical formulation of the equations, which are purely convective, and their numerical resolution, based on a finite-volume integration scheme. Even if the equations are not diffusive, the numerics generates its own diffusion, this contribution radically altering the model predictions. Diffusive fluxes in the population balance equation yield source terms in the DQMOM transport equations that, if neglected, lead to wrong results; therefore, for the model to work, we need to account for these terms. This issue cannot be circumvented; however, for populations where physical processes change the particle size continuously and discontinuously, the source terms related to these processes are likely to outweigh those related to numerical diffusion; if this is the case, the latter should not affect the numerical predictions significantly. The analysis presented in this work shows that purely convective quadrature transport equations cannot predict micromixing, even if each quadrature class is convected with its own velocity field. This problem concerns primarily the quadrature nodes, which result to be passive, nondiffusive scalars convected with the same velocity field as the quadrature weights. At least their transport equations should feature diffusion. Nevertheless, we cannot justify theoretically why such terms should arise. When modeling multiphase (or multicomponent, single-phase) systems we can choose between two modeling approaches: (1) solving for each phase (or component) a continuity equation and a dynamical equation, considering for each phase a different velocity field, or (2) solving for each phase a continuity equation and, assuming that all the phases share the same velocity, solving just one dynamical equation for the mixture. Only in the second approach the continuity equations feature diffusive fluxes.53 Now, since we have followed the first approach, diffusion should not arise. This problem, which has revealed itself because of the simplicity of the system examined, clearly needs to be researched further. Literature Cited (1) Sinclair, J. L.; Jackson, R. Gas-particle flow in a vertical pipe with particle-particle interactions. AIChE J. 1989, 35, 1473. (2) Ding, J.; Gidaspow, D. A bubbling fluidization model using kinetic theory of granular flow. AIChE J. 1990, 36, 523. (3) Kuipers, J.; van Duin, K.; van Beckum, F.; van Swaaij, W. Computer simulation of the hydrodynamics of a two-dimensional gas-fluidized bed. Comput. Chem. Eng. 1993, 17, 839. (4) Hrenya, C. M.; Sinclair, J. L. Effects of particle-phase turbulence in gas-solid flows. AIChE J. 1997, 43, 853. (5) Pain, C. C.; Mansoorzadeh, S.; de Oliveira, C. R. E. A study of bubbling and slugging fluidised beds using the two-fluid granular temperature model. Int. J. Multiphase Flow 2001, 27, 527. (6) Lettieri, P.; Cammarata, L.; Michale, G.; Yates, J. G. CFD simulations of gas-fluidized beds using alternative Eulerian-Eulerian modelling approaches. Int. J. Chem. React. Eng. 2003, 1, A5.

5152


(7) Cammarata, L.; Lettieri, P.; Michale, G.; Colman, D. 2D and 3D CFD simulations of bubbling fluidized beds using Eulerian-Eulerian models. Int. J. Chem. React. Eng. 2003, 1, A48. (8) Gelderbloom, S. J.; Gidaspow, D.; Lyczkowski, R. W. CFD simulations of bubbling/collapsing fluidized beds for three Geldart groups. AIChE J. 2003, 49, 844. (9) Lettieri, P.; Saccone, G.; Cammarata, L. Predicting the transition from bubbling to slugging fluidization using computational fluid dynamics. Chem. Eng. Res. Des. 2004, 82, 939. (10) Mazzei, L.; Lettieri, P. CFD simulations of expanding/contracting homogeneous fluidized beds and their transition to bubbling. Chem. Eng. Sci. 2008, 63, 5831. (11) van Wachem, B. G. M.; Schouten, J. C.; van den Bleek, C. M.; Krishna, R.; Sinclair, J. L. CFD modeling of gas-fluidized beds with a bimodal particle mixture. AIChE J. 2001, 47, 1292. (12) Wirsum, M.; Fett, F.; Iwanowa, N.; Lukjanow, G. Particle mixing in bubbling fluidized beds of binary particle systems. Powder Technol. 2001, 120, 63. (13) Howley, M. A.; Glasser, B. J. Hydrodynamics of a uniform liquidfluidized bed containing a binary mixture of particles. Chem. Eng. Sci. 2002, 57, 4209. (14) Huilin, L.; Yurong, H.; Gidaspow, D.; Lidan, Y.; Yukun, Q. Size segregation of binary mixture of solids in bubbling fluidized beds. Powder Technol. 2003, 134, 86. (15) Gera, D.; Syamlal, M.; O’Brien, T. J. Hydrodynamics of particle segregation in fluidized beds. Int. J. Multiphase Flow 2004, 30, 419. (16) Cooper, S.; Coronella, C. J. CFD simulations of particle mixing in a binary fluidized bed. Powder Technol. 2005, 151, 27. (17) Sun, Q.; Lu, H.; Liu, W.; He, Y.; Yang, L.; Gidaspow, D. Simulation and experiment of segregating/mixing of rice husk-sand mixture in a bubbling fluidized bed. Fuel 2005, 84, 1739. (18) Huilin, L.; Yunhua, Z.; Ding, J.; Gidaspow, D.; Wei, L. Investigation of mixing/segregation of mixture particles in gas-solid fluidized beds. Chem. Eng. Sci. 2007, 62, 301. (19) Owoyemi, O.; Mazzei, L.; Lettieri, P. CFD modeling of binaryfluidized suspensions and investigation of role of particle-particle drag on mixing and segregation. AIChE J. 2007, 53, 1924. (20) Mathiesen, V.; Solberg, T.; Hjertager, B. H. An experimental and computational study of multiphase flow behavior in a circulating fluidized bed. Int. J. Multiphase Flow 2000, 26, 387. (21) Mathiesen, V.; Solberg, T.; Hjertager, B. H. Predictions of gas/ particle flow with an Eulerian model including a realistic particle size distribution. Powder Technol. 2000, 112, 34. (22) Syamlal, M.; Rogers, W. A.; O’Brien, T. J. MFIX Documentation and Theory Guide; DOE/METC94/1004, NTIS/DE94000087; 1993; electronically available from http://www.mfix.org. (23) Gidaspow, D. Multiphase Flow and Fluidization; Academic Press: New York, 1994. (24) Drew, D. A.; Passman, S. L. Theory of Multicomponent Fluids; Applied Mathematical Sciences; Springer: New York, 1998. (25) Olmos, E.; Gentric, C.; Vial, C.; Wild, G.; Midoux, N. Numerical simulation of multiphase flow in bubble column reactors. Influence of bubble coalescence and break-up. Chem. Eng. Sci. 2001, 56, 6359. (26) Lehr, F.; Mewes, D. A transport equation for the interfacial area density applied to bubble columns. Chem. Eng. Sci. 2001, 56, 1159. (27) Buwa, V. V.; Ranade, V. V. Dynamics of gas-liquid flow in a rectangular bubble column: experiments and single/multi-group CFD simulation. Chem. Eng. Sci. 2002, 57, 4715. (28) Venneker, B. C. H.; Derksen, J. J.; van den Akker, H. E. A. Population balance modeling of aerated stirred vessels based on CFD. AIChE J. 2002, 48, 673. (29) 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.

(30) Fan, R.; Fox, R. O. Segregation in polydisperse fluidized beds: Validation of a multi-fluid model. Chem. Eng. Sci. 2008, 63, 272. (31) Petrovsky, I. Lectures on the Theory of Integral Equations; Graylock Press: Rochester, 1957. (32) Ramkrishna, D. Population Balances; Academic Press: New York, 2000. (33) Zucca, A.; Marchisio, D. L.; Vanni, M.; Barresi, A. A. Validation of bivariate DQMOM for nanoparticle processes simulation. AIChE J. 2007, 53, 918. (34) Randolph, A. D.; Larson, M. A. Theory of particulate processes; Academic Press: New York, 1971. (35) Hulburt, H. M.; Katz, S. Some problems in particle technology. Chem. Eng. Sci. 1964, 19, 555. (36) Marchisio, D. L.; Fox, R. O. Solution of population balance equations using the direct quadrature method of moments. J. Aerosol Sci. 2005, 36, 43. (37) McGraw, R. Description of aerosol dynamics by the quadrature method of moments. Aerosol Sci. Technol. 1997, 27, 255. (38) Gordon, R. G. Error bounds in equilibrium statistical mechanics. J. Math. Phys. 1968, 9, 655. (39) Dette, H.; Studden, W. J. The theory of canonical moments with applications in statistics, probability, and analysis; Wiley: New York, 1997. (40) Geldart, D. Gas Fluidization Technology; Wiley: New York, 1973. (41) Marchisio, D. L.; Fox, R. O. Multiphase reacting flows: modelling and simulation; CISM Courses and Lectures 492; International Centre for Mechanical Sciences; Springer: New York, 2007. (42) Lettieri, P.; Mazzei, L. Challenges and issues on the CFD modeling of fluidized beds: a review. J. Comput. Multiphase Flows 2009, 1, 83. (43) Jackson, R. The Dynamics of Fluidized Particles; Cambridge Monographs on Mechanics; Cambridge University Press: New York, 2000. (44) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena; Wiley: New York, 2007. (45) Schaeffer, D. G. Instability in the evolution equations describing incompressible granular flow. J. Differ. Equations 1987, 66, 19. (46) Mazzei, L.; Lettieri, P. A drag force closure for uniformly-dispersed fluidized suspensions. Chem. Eng. Sci. 2007, 62, 6129. (47) Dallavalle, J. M. Micromeritics; Pitman: New York, 1948. (48) Richardson, J. F.; Zaki, W. N. Sedimentation and fluidization: Part I. Trans. Inst. Chem. Eng. 1954, 32, 35. (49) Rowe, P. N. A convenient empirical equation for estimation of the Richardson & Zaki exponent. Chem. Eng. Sci. 1987, 42, 2795. (50) Syamlal, M. The particle-particle drag term in a multiparticle model of fluidization; National Technical Information Service; DOE/MC/213532373, NTIS/DE87006500; 1987. (51) Lo, S. M. Mathematical basis of a multiphase flow model. U.K. At. Energy Auth., Harwell Lab., [Rep] AERE-R 1989, AERE-R13432. (52) Mazzei, L. Eulerian modelling and computational fluid dynamics simulation of mono and polydisperse fluidized suspensions. Ph.D. Dissertation, Department of Chemical Engineering, University College London, 2008. (53) Curtiss, C. F.; Bird, R. B. Multicomponent diffusion in polymeric liquids. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 7440. (54) Fox, R. O. Computational Models for Turbulent Reacting Flows; Cambridge University Press: New York, 2003. (55) Ferziger, J. H.; Peric, M. Computational methods for fluid dynamics; Springer: New York, 2002.

ReceiVed for reView July 10, 2009 ReVised manuscript receiVed November 6, 2009 Accepted November 13, 2009 IE901116Y

Direct Quadrature Method of Moments for the Mixing of Inert

Recommend Documents