Prediction of Radial Distribution Function of Particles in a Gas− Solid

Ind. Eng. Chem. Res. 2009, 48, 1343–1352

1343

Prediction of Radial Distribution Function of Particles in a Gas-Solid Fluidized Bed Using Discrete Hard-Sphere Model Shuyan Wang,† Guodong Liu,† Huilin Lu,*,† Bai Yinghua,† Jianmin Ding,‡ and Yunha Zhao† School of Energy Science and Engineering, Harbin Institute of Technology, Harbin, 150001, China, and Agere Systems Inc., 555 Union BouleVard, Allentown, PennsylVania 18109-3286

Flow behavior of particles in a two-dimensional bubbling fluidized bed is predicted by using discrete hardsphere model for particle-particle collision. Quantities of radial distribution functions of monosized particles, binary-sized particles, and binary density particles are obtained. Experimentally or theoretically proposed formulations for the radial distribution functions are evaluated based on our numerically predicted results. For monosized particles, the simulated radial distribution functions are in agreement with computed results from both Bagnold equation (1954) and Ma and Ahmadi (1986) equation. For the binary mixture with the different diameters but identical density, the pair radial distribution functions proposed by Boublik (1970) and Mansoori et al. (1971) agree with simulated data. For binary mixture of different densities, a modified equation of the pair radial distribution function is proposed to correlate from our simulation results. 1. Introduction The kinetic theory of granular flow as originated in the middle of 19801-3 has been developed and widely used for numerical simulations of gas-solid fluidization since early 1990.4-11 In the kinetic theory model, the transportation equations of solid phase in the gas-solid two-phase flow are derived from the Boltzmann equation with the binary particle collision term expressed by a pair distribution function. From Enskog assumption, the pair distribution function is proportional to a product of two particles distribution function involved in the collision. Based on the kinetic theory of dense gases,12 the product factor is the radial distribution function at the point of two particles’ contact. In general, the radial distribution function denoted as go is a measure of probability of binary particle collisions and can be a function of particle concentration.13 The quantity of the radial distribution function is equal to unity as the particle concentration reduces to zero and increases with the increase of particle concentration, and becomes infinity as particles are closely packed. Constitutive equations of solid phase in gas-solid flow such as solid viscosities and fluctuation energy dissipations derived from kinetic theory of granular flow are directly related to the radial distribution function. Therefore, results from numerical simulation of gas-solid flow in fluidized bed using the kinetic theory of granular flow model are affected by the values of radial distribution function. Radial distribution functions play a key role in the statistical mechanics of liquid.14 They can be measured using fluorescence microscopy15 and neutron scattering techniques.16 The radial distribution function of different liquids is computed by Monte Carlo (MC) molecular simulations.17,18 The radial distribution function is computed from molecular dynamics (MD) simulations to examine the nonequilibrium structure formation in a colloidal suspension.19 The radial distribution function of fluid catalytic cracking (FCC) particles is calculated from the measured local particle distributions by means of a chargecoupled device (CCD) camera system in a circulating fluidized bed.20 * To whom correspondence should be addressed. Tel.: +0451 8641 2258. Fax: +0451 8622 1048. E-mail: [email protected]. † Harbin Institute of Technology. ‡ Agere Systems Inc.

In present study, we use the discrete particle model to compute the particle motion. The particle-particle collisions and external forces acting on the particles are directly taken into account. In general, the mechanism of particle-particle collisions can be described by soft-sphere or hard-sphere models.21 In the soft-sphere model, the Newtonian equation of motion of individual particles is integrated. The collisions between particles are simulated using Hooke’s linear springs and dash pots. In this approach, the particles are allowed to be slightly overlapped and the contact forces are subsequently calculated. This model has been previously used for investigation of interparticle force effect on fluidization characteristics in the fluidized beds.22-25 On the other hand, the hard-sphere models have been used to model particle-particle collisions.26-28 The interactions between particles are assumed to be pairwise additive and instantaneous. In the simulation, the collisions are processed one by one according to the order in which the events occur. The motion of particles can be predicted in the fluidized bed.29,30 Note that the possible occurrence of multiple collisions at the same instant cannot be accounted for. In present investigation, we use the hard-sphere discrete particle model to compute flow behavior of monosized particles and polydispersed particles with different diameters or densities in a bubbling fluidized bed. The radial distribution functions of these types of particles are predicted. Comparisons between the present results and the theoretically proposed formulations used in the literature are made. The equation of radial distribution function at contact with the different densities of particles is proposed. 2. Eulerian-Lagrangian Gas-solid Flow Model 2.1. Continuity and Momentum Equation for Gas Phase. The Eulerian-Lagrangian method solves simultaneously the Navier-Stokes equation for the gas phase and the Newtonian equations for particle motions. For the gas phase, the equations of mass and momentum are21 ∂ (ε F ) + ∇ ·(εgFgvg) ) 0 ∂t g g

10.1021/ie8007049 CCC: $40.75  2009 American Chemical Society Published on Web 12/29/2008

(1)

1344 Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009

∂ (ε F v ) + ∇ ·(εgFgvgvg) ) -εg ∇ p + ∇ ·τg + εgFgg ∂t g g g

Np(V)

∑f

i

i)1

(2) where g is the acceleration due to gravity, p the thermodynamic pressure, εg porosity, Fg gas density, and τg the viscous stress Np(V) tensor. The coupling term ∑i)1 fi between gas phase and solid phase is estimated as the sum of the drag on each particle within the corresponding fluid control volume. The stress tensor of gas phase can be represented as 2 (3) τg ) µg[∇vg + (∇vg)T] - µg(∇·vg)I 3 where µg is the viscosity of gas phase. 2.2. Equations for Particles. Analysis of particles collision in this study is based on collision dynamics with the following assumptions: (1) particles are spherical and quasi-rigid; (2) collisions are binary and instantaneous; (3) interaction forces are impulsive, and all other infinite forces are negligible during collision;. (4) particle motion is two-dimensional with the particle mass center moving in one plane; (5) particles are in free flight in between collisions; and (6) both the coefficient of restitution and the friction coefficient are constant in the simulations. For binary collision of particles, the following equations can be derived by applying Newton’s second and third laws:21,26,27 m1(u1 - u1,0) ) J

(4)

m2(u2 - u2,0) ) -J

(5)

2

m1d1 (ω1 - ω1,0) ) n × J 4

(6)

m2d22 (ω2 - ω2,0) ) -n × J (7) 4 where d and m are particle diameter and mass, respectively. u and ω are respectively translational velocity and angular velocity of particles. J is the impulse force. From eqs 4-7 it is clear that the postcollision velocities of both particles can be calculated when the impulse is known. The relative velocity u12 at the point of contact is u12 ) u1 - u2

(8)

The tangential unit vector is t)

u12,0 - n(u12,0·n) |u12,0 - n(u12,0·n)|

contact and the collision is of the sliding type. The nonsliding collisions are of the sticking type. When the coefficient of tangential restitution is equal to zero, the tangential component of the relative velocity becomes zero during a sticking collision. When the coefficient of tangential restitution is greater than zero in such a collision, reversal of the tangential component of the relative velocity will occur. The criterion to determine the type of collision on basis of precollision information is as follows:21,26 µ
0 these particles are moving away from each other and they will not collide. For the case of a collision between particle and wall, the collision time depends simply on the distance of the particle to the wall and on the normal velocity component toward the wall. However, this equation does not consider the effect of particle acceleration during the collisions. For the binary mixture system with different diameter or density of particles, the calculated collision time between two particles includes the different accelerations of particles due to the different diameter or density. Considering two spherical particles located at r1 and r2 whose velocities are W1 and W2 experiencing different accelerations a1 and a2, respectively, these particles will collide at time t + t12, namely |r12 + V12t12 + (1/ 2 | ) (1/2)(d1 + d2). Then, t12 may be obtained by finding 2)a12t12 the smallest real positive root of the following quadratic equation:34,35

Figure 1. Adjacent PCVs included in the collision. Table 1. Definition of Variables Used in the Solution of the Quadratic Equation D1 ) (V212 + δ12)2 - 6χ12ζ12 + 3a212[r212 - 0.25(d1 + d2)2] D2 ) 2(V212 + δ12)3 -18χ12ζ12(V212 + δ12) + 27a212ζ212 + 9 [r212 - 0.25(d1 + d2)2][3χ212 - 2a212(V212 + δ12)] D3 ) [-4D31 + D22]0.5 D4 ) 4[(χ12/a12)2 - (2/3)(V212 + δ12)]/a212 D5 ) 25/3[(22/3D1)/(D2 + D3)1/3 + (D2 + D3)1/3]/3a212 D6 ) 16[χ12(V212 + δ12)/a212 - 4(ζ12 + χ312/a412)]/a212

(T1-1) (T1-2) (T1-3) (T1-4) (T1-5) (T1-6)

1 2 4 1 a t + χ12t312 + (V212 + δ12)t212 + 2ζ12t12 + r212 - (d1 + d2)2 ) 4 12 12 4 0 (24) Here, a12 and V12 are the values of the relative acceleration and velocity vectors between particles 1 and 2, respectively. The parameter χ12 is the dot product of these two vectors χ12 ) a12 · W12, and the variable r12 denotes the distance between the centers of the particles. The parameter δ12 ) a12 · r12 represents the dot product of relative acceleration and position vectors. The dot product of the relative velocity and relative position vectors of particle 1 and 2 is ζ12, i.e., ζ12 ) r12 · W12. The solutions of eq 24 can be explicitly expressed by

[

t12 ) -

[

t12 ) -

χ12 a212

χ12 a212

] [

1 1 - (D4 + D5)0.5 - 2D4 - D5 2 2

] [

Figure 2. Principle of area-weighted averaging. Table 2. Parameters Used in Simulations of Monosized Particles

]

-0.5 0.5

D6(D4 + D5) 4

1 1 + (D4 + D5)0.5 - 2D4 - D5 + 2 2 D6(D4 + D5)-0.5 4

]

(25a)

0.5

(25b)

where the expressions for the variables D4, D5, and D6 are listed in Table 1. It is noted that neither one nor both pairs of roots may be real. The smallest positive real root of the quadratic eq 25a corresponds to impact. Obviously, the absence of any real positive root implies that the pair of particles will not imminently collide. 2.5. Neighbor List. In molecular dynamic simulations, the nearest-neighbor list is constructed by defining a cutoff distance Rlist which gives a circular area in the two-dimensional system (or a spherical space in the three-dimensional system). This involves the computation of circulars or spheres at the cost of CPU time. Hoomans et al.26 used the square neighbor list area in the two-dimensional system. All particles scanned within the square of size Dlist are listed. The size of Dlist is chosen to be 5 times the particle diameter. In this study, we have developed an algorithm for detecting particle binary collision in order to save CPU time. Several fluid control volumes (FCV) are taken to represent an elementary volume in the Eulerian coordinate

diameter, d (mm) density, Fs (kg/m3) coeff of restitution, e dynamic friction coeff, µ coeff of tangential restitution no. of particles, n

2.3 2600 0.95 0.1 0.3 1923

width, B (mm) height, H (mm) cell width, ∆x (mm) cell height, ∆y (mm) no. x-cell no. y-cell

150 650 10 20 15 25

system, shown in Figure 1. Each particle in the FCV represents a particle control volume (PCV). Every particle within the FCV (particle i) and its adjacent four PCVs are included in a list, so that this collision searching is only limited to the particle i located FCV and its nearest four neighbor FCVs (the hatched parts in Figure 1). We try to include as many PCVs within an FCV as possible in order to make the collision list as small as possible. For the first PCV in the collision list, we search for overlapping of the particle with other particles in its own FCV and the nearest four neighbor FCVs. When a first overlapping is detected, we treat the collision dynamics for the two particles involved. For the first particle the search for overlapping is terminated. Then, we repeat the searching procedure for the next particles in the list. For these particles, we only search for overlapping with the subsequent particles. When all of PCV within the respective FCV have undergone the procedure, we repeat the same treatment for the next PCV until every collision list for all the PCVs in the calculation domain has been treated. Nevertheless, this method seems to be efficient as no insignificant fraction of particle overlapping is observed even in high particle concentration zones.

1346 Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009

2.6. Grid Mapping. For the calculation of the force acting on a suspended particle, the local averaged values of pressure, porosity, and velocities of gas phase at the position of the particle (the Lagrangian grid) are required. From the numerical solution method used, these variables are only known at discrete nodes of the computational domain. An area-weighted averaging j of a technique is used to obtain the local averaged value Q quantity Q(i,j) from the four surrounding computational nodes as shown in Figure 2. The local averaged value can be calculated as follows:26-28,36 j EfL ) Q

Ai,jQi,j + Ai+1,jQi+1,j + Ai,j+1Qi,j+1 + Ai+1,j+1Qi+1.j+1 ∆x∆y (26) Figure 4. A close look of a part of instantaneous particle distribution in the fluidized bed.

with Ai,j ) δxδy Ai+1,j ) (∆x - δx)δy Ai,j+1 ) δx(∆y - δy) Ai+1,j+1 ) (∆x - δx)(∆Y - δy) The distances δx and δy in this averaging technique are calculated from the position of the particle in the staggered grid. 2.7. Calculation of Porosity. For each cell of the computational domain, porosity can be calculated on the basis of the area occupied by the particles in the cell. The porosity is obtained by subtracting the sum of the particle volumes from the volume of a fluid cell. When a particle overlaps one neighboring cell or more, the volume fraction included in a cell is taken into account to calculate the porosity of the cell. A “2D” porosity in the cell is defined as n

∑s

Figure 5. Calculation of radial distribution function.

i

εg,2D ) 1 -

i)1

∆x∆y

(27)

where ∆x and ∆y are the widths of computational cell and s is the surface area of the particle located in the cell. However, the above-defined “2D” porosity is inconsistent with the applied empiricism in the calculation of the drag force exerted on a particle and of the interfacial friction. Since the relevant correlations are from 3D systems. Thus, we used a 3D porosity εg,3D26-28,36 εg,3D ) 1 -

2

√π√3

(1 - εg,2D)3⁄2

(28)

2.8. Simulation Procedures. The particle motion is divided into two steps: collision step and free flight step. The former is

controlled by dynamics of particle collision, and the latter is controlled by dynamics of particle motion. The simulation procedures are (1) to randomly set particles in the bed to obtain initial condition, then gas flow is uniformly induced into the bed; (2) to calculate drag force acting on a particle; from the position and velocity of particles, a sequence of collisions is processed using the collision model; at the end of this step, the new position and velocity of particles are determined; the porosity and force acting on the gas of a cell are obtained; and (3) to obtain gas flow field by solving eqs 1 and 2. In all of the following simulations, the computer time step for gas phase is set to be 10-5 s. The total time of simulations are continued for 10 s.

Figure 3. Simulations of bubble and particles motions in a fluidized bed with coefficient of restitution of 0.95 at the superficial gas velocity of 2.5 m/s.

Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 1347

Figure 6. Calculated radial distribution function of monosized particles in bubbling fluidized bed.

through the bed and breaks at the bed surface. Such similar simulations are also found in the published literature.26-30 3.2. Radial Distribution Function of Monosized Particles. In molecular dynamics simulations, the radial distribution function is often applied to characterize the structure of atomic fluids. In the kinetic theory of granular flow, as discussed earlier, the radial distribution function is a function of the volume occupied by the particles. The configurational distribution function defines a FN(r1,r2,...rN), which is a multidimensional probability density function.38,39 It represents the probability per unit volumeN finding N particles at locations (r1,r2,...rN) where ri represents the position vector of particle i. The integral of FN(r1,r2,...rN) over all possible locations of these N particles yields unity

∫ F (r , r , r ...r ) dV V

Figure 7. Calculated instantaneous radial distribution function at contact and concentration of particles.

N

1

2

3

N

1

dV2 dV3...dVN ) 1

(29)

If the integration is performed over all but one set of particle coordinates, ∫V FN(r1,r2,...rN) dV2...dVN, the resulting integral represents the probability per unit volume of finding a particle at location r1. The single-particle distribution function, F1(r1), is F1(r1) ) N

∫ F (r , r , r ...r ) dV V

N

1

2

3

N

2

dV3...dVN ) n(r1) (30)

where n(r1) is the local number density (the number of particles per unit volume at r1). The two-particle distribution function is obtained by integrating FN(r1,r2,...rN) over all but two sets of particle coordinates F2(r1, r2) ) N(N - 1)

∫ F (r , r , ...r ) dV ...dV V

N

1

2

N

3

N

(31)

The radial distribution function g(r) is defined as g(r) )

Figure 8. Calculated radial distribution function of monosized particles as a function of concentrations.

3. Results and Discussion 3.1. Fluidization of Monosized Particles. The parameters used for fluidization of monosized particles are listed in Table 2. The minimum fluidization velocity of particles is 1.77 m/s based on the equation proposed by Wen and Yu.37 The calculated terminal velocity of particle is 16.5 m/s. The maximum particle volume fraction at maximum packing conditions is assumed to be 0.64. Snapshots of the simulated particle motion are shown in Figure 3 at the superficial gas velocity of 2.5 m/s. Particles are fluidized by gas flow in the bed. Simulation results show that a bubble is initiated at the orifice and gradually grows in size

F2(r1, r2)

(32) n2 where r ) r2 - r1 is the separation vector between the two particles and n is the average number density. Figure 4 shows a typical computed spatio-temporal particle distribution. Coordinates of each particle are recorded. From the spatial distribution of particles, the radial distribution function g(r) can be calculated on the basis of eq 32: ∆N × AREA (33) 2πr∆rN where N is the total particle number in the AREA, and ∆N is the particle number in the computing area.38,39 Figure 5 gives the concept of the radial distribution function calculation. As r f ∞, local density ) local density × g(r). Hence, g(r) f 1.0. The typical radial distribution function profiles in the fluidized bed are shown in Figure 6. The highest quantity of radial g(r) )

1348 Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009

Figure 9. Simulations of bubble and binary mixture with the small (red) and large (blue) particle diameters in a fluidized bed at the superficial gas velocity of 2.5 m/s.

Figure 10. Calculated pair radial distribution function of polydispersed particles with different diameters and the same densities in bubbling fluidized bed. Table 3. Parameters Used in Simulations of Polydispersed Particles with the Different Diameters diameter of smaller particle, d1 (mm) diameter of larger particle, d2 (mm) density, Fs (kg/m3) no. of larger particles, n2 coeff of restitution, e dynamic friction coeff, µ coeff of tangential restitution

2.3 4.0 2600 500 0.95 0.1 0.3

distribution function is at the distance r of a particle diameter, then decreases to unity as the dimensionless distance r/d increases. The radial distribution function is a local measure of how close the particle distribution is to uniform (a distribution where all particle locations are equally likely). g(r) ) 1.0 implies that the probability of finding a particular particle in the subvolume dV at r is the same as would be obtained if the particles were uniformly distributed. As g(r)is greater or less than unity, it implies an enhanced or reduced probability relative to the uniform distribution. From these figures, the radial distribution function at contact (r ) d), go, is determined. Figure 7 shows the distributions of radial distribution function at contact and concentration of particles as a function of time. The local concentration of particles and the radial distribution function are oscillations with time. The high concentration of particles gives a high radial distribution function at contact. Thus, the radial distribution function could be mainly related with particle concentrations. Figure 8 shows the distributions of radial distribution function at contact as a function of concentration of particles. The radial distribution function increases with the increase of particle concentrations. The radial distribution function at contact was proposed by Carnahan and Starling40 for dense rigid spherical gases: go(εs) )

2 - εs 2(1 - εs)3

(34)

width, B (mm) height, H (mm) no. of smaller particles, n1 cell width, ∆x (mm) cell height, ∆y (mm) no. x-cell no. y-cell

150 650 1000 10 20 15 25

where εs is concentration of particles. Values from above equation are close to the results from molecular dynamics simulations for concentration of particles up to 0.55. Savage41 proposed a simple expression that was implicit in the work of Bagnold:42

[ ( )]

go(εs) ) 1 -

εs

1⁄3 -1

εs,max

(35)

Savage41 noticed that the predictions of stresses were close to experimental results if the maximum packing density for randomly packed spheres, εs,max ) 0.64, was applied. From the kinetic theory of dense gases, Ma and Ahmadi43 proposed a formula of radial distribution function at contact for hard particles go(εs) ) 1 + 4εs

1 + 2.5εs + 4.5904εs2 + 4.515439εs3 (36) εs 3 0.67802 1εs,max

[ ( )]

The formula results from eqs 34, 35, and 36 are shown in Figure 8. We see that the calculated radial distribution function by Carnahan and Starling equation is lower than that by present simulations. The quantity of the radial distribution function from eq 34 does not go to infinite as particle concentrations approach to the maximum packing value. We see that the present predicted

Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 1349

Figure 11. Profile of radial distribution function of small and large particles in binary mixture with different diameters and same density of particles in bubbling fluidized bed.

Figure 12. Simulations of bubble and binary mixture with the light (red) and heavy (blue) particle densities in a fluidized bed at the superficial gas velocity of 2.5 m/s.

results are in agreement with those from eqs 35 and 36 in the bubbling fluidized bed of monosized particles. 3.3. Radial Distribution Function of Binary-Sized Particles. Flow behavior of binary-sized particles with the same density is numerically predicted in the fluidized bed. In a certain velocity range, segregation can be displayed with small particles (flotsam) in the top part and larger particles (jetsam) at the bottom part of the bed. The parameters used are listed in Table 3. The diameters of particles are 2.3 and 4.0 mm with the same density of 2600 kg/m3. Figure 9 shows the instantaneous solid flow patterns at the superficial gas velocity of 2.5 m/s. The force of gas acting on a particle is varied for different sized particles. Because the drag force is inversely proportional to the particle diameter, more small particles than large particles are carried up by the incoming gas. This leads to binary mixture segregation by fluidizing gas. Here we define pair radial distribution functions, one is the self-pair radial distribution function, gii, and the other is the mutual radial distribution function gij for the radial distribution between i type of particles and j type of particles.44 During the

computations, the distances between the two types of particles in the cells are counted in the same way as calculations in the monocomponent particles system. Figure 10 shows the pair radial distributions of binary mixture with the different diameters at the superficial gas velocity of 2.5 m/s. The quantities of the radial distribution function decrease with the increase of dimensionless distance. The first peak occurs at a radius close to the size parameter for each component in the mixture (g11(r) and g22(r)). The peak of g12(r) occurs at the dimensionless distance of r/((d1 + d2)/2) ) 1.0. Figure 11 shows the self-pair and mutual pair radial distribution functions at contact as a function of concentrations. Similarly, all quantities of these pair radial distribution functions increase with the increase of particle concentrations. An extension of the equations of state of Carnahan and Starling to hard-sphere mixture may be obtained and leads to a so-called Boublik-Mansoori-Carnahan-Starling-Leland equation of state. The pair radial distribution function proposed by Boublik45 and Mansoori et al.46 takes the form

1350 Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009

go,Rβ )

(

)

dRdβ η2 1.5η3 dRdβ η2 2 0.5η32 1 + + 1 - η3 dRβ η3 (1 - η )2 dRβ η3 (1 - η )3 3 3 (37)

where dRβ ) (dR + dβ)/2, and η is the generalized packing factor defined as 2

ηk )

∑ π6 j)1

εs,jFs,j k d (k ) 2, 3) ms,j j

[ ( )]

go,Rβ ) 1 -

and εs,j is the concentration of particles species j. This equation was used by Huilin et al.47 to calculate pair radial distribution function in a bubbling fluidized bed of binary particles mixture. On the basis of the single solid-phase model of eq 35, Mathiesen et al.48 proposed a binary radial distribution function in the numerical simulations of fluidized beds:

1⁄3 -1 (ε + ε ) sR sβ

1 - εg

(39)

Recently, Fan and Fox49 used the following equation to predict the pair radial distribution function at contact over two components in the bubbling fluidized beds: go,Rβ )

(38)

1 - εg εs,max

2 3dRdβ εsλ 1 + 2 εg (d + d )ε λ)1 dλ R β g

∑

(40)

2 where εg is porosity, εg ) 1 - ∑λ)1 εs,λ. This equation can be directly deduced from eq 37 as the third term of right-hand side of eq 37 is neglected. The calculated pair radial distributions from eqs 39 and 40 are also shown in Figure 11. The pair radial distribution functions from the present simulations and calculated from the above equations give similar increments with the particle concentrations. Data from present simulations and

Figure 13. Profile of radial distribution function of light and heavy particles in the binary mixture with different densities and same diameter of particles in bubbling fluidized bed.

Figure 14. Profile of radial distribution function of mixture particles with different densities and same diameter of particles in bubbling fluidized bed.

Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 1351

from eqs 37 and 40 are close to each other. But results from eq 39 show much higher values at high particle concentrations. It is also noted that there are some unrealistic values. The calculated self-pair radial distribution function go,11 from eq 39 is less than unity at the low concentrations, since the viscosity and collisional energy dissipation of the solid phase are proportional to the quantity of radial distribution function. Thus, the radial distribution functions are directly affecting the determination of energy exchange and dissipations of solid phase in the numerical simulations of fluidized beds. 3.4. Radial Distribution Function of Particles with Two Different Densities. Fluidization of particles with two different densities and same diameter of particles in the bubbling fluidized bed is numerically predicted. The bed consists of 960 light (flotsam) particles with a density of 1200 kg/m3 and 960 heavy (jetsam) particles with a density of 2600 kg/m3. These particles have an identical diameter of 2.3 mm. The restitution coefficient, dynamic friction coefficient, and coefficient of tangential restitution of both light and heavy particles are 0.9, 0.1, and 0.3, respectively. The superficial gas velocity is 2.5 m/s. Figure 12 shows the instantaneous particles distributions at several times in the bed. The driving force for particle motions is the gas fluidization. Initially, particles are well mixed and settled in the bed. The upward injecting gas drags the particles in the gas flow direction to the top of bed by overcoming the particle gravity. Even under the same drag force, those lighter particles can be easier carried by gas to the top of the bed, and leading to more light particles floated at the top of the bed, while more heavy particles found near the bottom of the bed. Such phenomena are shown in Figure 12. Figure 13 shows the pair radial distribution functions of the two different particle densities at the superficial gas velocity of 2.5 m/s in bubbling fluidized bed. Although results from the present simulations and from eqs 37, 39, and 40 increase with the particles’ concentrations, eq 39 gives some unreasonable values of the self-pair radial distributions of light and heavy particles at the low concentrations, that is, similar to the case of binary particles presented in the previous section. It is also noted that the self-pair radial distribution functions at contact, go,11 and go,22 from eqs 37 and 40, are in agreement with present simulations, while the calculated mutual pair radial distribution functions go,ij from eqs 37 and 40 are higher than that present simulations. Considering the different densities of particles in the binary mixture, we modify eq 37 and propose following equation for pair radial distribution function at contact go,Rβ )

(

1 3 dRdβ + 2 εg ε dR + dβ g

)∑

(

εsλ(FsR + Fsβ) 0.5 dRdβ + 3 2Fsλdλ ε dR + dβ λ)1 2

(∑

g

)

2

)

εsλ(FsR + Fsβ) 2 (41) 2Fsλdλ λ)1 2

The mutual pair radial distribution functions at contact calculated from eq 41 are given in Figure 13. It can be seen that the calculated mutual pair radial distribution functions at contact from eq 41 are in agreement with the present simulations. Figure 14 shows the pair radial distribution functions of the two different particle densities with the same particle diameter of 1.0 mm at the superficial gas velocity of 2.0 m/s in a bubbling fluidized bed. The pair radial distribution functions increase with the increase of particle concentrations. At the low concentrations of particles, eq 39 gives some unrealistic values of the self-pair radial distributions of light and heavy particles. The predicted mutual pair radial distribution functions at contact from eq 41 agree with present simulations.

4. Conclusions Accurate characterization of radial distribution function is important because constitutive equations of solid phase derived from kinetic theory of granular flow depend on the radial distribution function. The computed results of gas-solid flow in fluidized bed are therefore affected by the radial distribution function. Our present simulations using the hard-sphere discrete model predict the quantities of the radial distributions of monosized particles and binary mixture of different diameters or densities. Based on our numerically predicted results, previous experimentally or theoretically proposed formulations for the radial distribution functions are evaluated. For monosized particles, the simulated radial distribution functions are in agreement with computed results from both Bagnold equation42 and Ma and Ahmadi43 equation. For the binary mixture with the different diameters but identical density, the pair radial distribution function from eq 37 proposed by Boublik45 and Mansoori et al.46 agrees with present simulated data. For binary mixture of different densities, a modified pair radial distribution function eq 41 is proposed to correlate from our simulation results. However, one should keep in mind the limitation of the discrete hard-sphere model. Just as the kinetic theory of granular flow, it assumes that collisions are binary and instantaneous. This makes the model an appropriate tool to predict the flow behavior for dilute particle suspensions, but not such realistic representation of highly inelastic particulate systems. In such systems, multiparticle and long-term contacts will occur at high particle concentrations. Simulations reveal the radial distribution function is required to correlate the observed dependence of the particle concentration and collision parameters. This is a good basis for the start of the determination of radial distribution function in the fluidized beds but it certainly is an approximation. More detailed descriptions of the particle collisions are required. Experimental measurement of radial distribution function, although arduous, is necessary in order to generate a sound basis for model validation and will be further pursued. Acknowledgment This work was supported by Natural Science Foundation of China through Grant No. 50776023 and National Key Project of Scientific and Technical Supporting Programs Funded by Ministry of Science and Technology of China (NO.2006BAA03B01), and by the Scientific Research Foundation of Harbin Institute of Technology through Grant No. HIT. 2003.34. Literature Cited (1) Savage, S. B.; Jeffrey, D. J. The stress tensor in a granular flow at high shear rates. J. Fluid Mech. 1981, 110, 255. (2) Jenkins, J. T.; Savage, S. B. A theory for the rapid flow of identical, smooth, nearly elastic, spherical particles. J. Fluid Mech. 1983, 130., 187. (3) Lun, C. K. K.; Savage, S. B.; Jeffrey, D. J.; Chepurniy, N. Kinetic theories for granular flow: Inelastic particles in Couette flow and slightly inelastic particles in s general flow field. J. Fluid Mech. 1984, 140, 223. (4) Ding, J.; Gidaspow, D. Bubbling fluidization model using kinetic theory of granular flow. AIChE J. 1990, 36, 523. (5) Nieuwland, J. J.; van sint Annaland, M.; Kuipers, J. A. M.; van Swaaij, W. P. M. Hydrodynamic modeling of gas/particle flows in riser reactor. AIChE J. 1996, 42, 1569. (6) Benyahia, S.; Arastoopour, H.; Knowlton, T. M.; Massah, H. Simulation of particles and gas flow behavior in the riser section of a circulating fluidized bed using the kinetic theory approach for the particulate phase. Powder Technol. 2000, 112, 24. (7) Wei, W.; Li, Y. Hydrodynamic simulation of fluidization by using a modified kinetic theory. Ind. Eng. Chem. Res. 2001, 40, 5066.

1352 Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 (8) Hadinoyo, K.; Jennifer, C. S. Effect of interstitial fluid on particleparticle interactions in kinetic theory approach of dilute turbulent fluidparticle flow. Ind. Eng. Chem. Res. 2004, 43, 3604. (9) Andrews, A. T.; Loezos, P. N.; Sundaresan, S. Coarse-grid simulation of gas-particle flows in vertical risers. Ind. Eng. Chem. Res. 2005, 44, 6022. (10) Sun, J.; Battaglia, F. Hydrodynamic modeling of particle rotation for segregation in bubbling gas-fluidized beds. Chem. Eng. Sci. 2006, 61, 1470. (11) Vaishali, S.; Roy, S.; Bhusarapum, S.; Al-Dahhan, M. H.; Dudukovic, M. P. Numerical simulation of gas-solid dynamics in a circulating fluidized-bed riser with Geldart group B particles. Ind. Eng. Chem. Res. 2007, 46, 8620. (12) Chapman, S.; Cowling, T. G. The mathematical theory of nonuniform gases, 3rd ed.;, Cambridge University Press: Cambridge, UK, 1970. (13) Gidaspow, D. Multiphase flow and fluidization: Continuum and kinetic theory descriptions; Academic Press: Boston, 1994. (14) Hansen, J. P.; McDonald, I. R. Theory of simple liquids; Academic Press: New York, 1991. (15) Yoshida, H.; Ito, K.; Ise, N. Microscopic observation and quasielastic light-scatting measurements of colloid crystals: Determination of the radial distribution function and structure factor for the two-state structure. J. Am. Chem. Soc. 1990, 112, 592. (16) Xu, W.; Mikolov, A. D.; Wasan, D. T. Role of depletion and surface-induced structural forces in bidisperse suspensions. AIChE J. 1997, 43, 3215. (17) Heyes, D. M. The liquid state: applications of molecular simulations; J. Wiley & Sons: Chichester, UK, 1997. (18) Wong, A. C. T.; Yu, K. W. Heterogeneous aggregation in binary colloidal alloys. Physica A 2002, 312, 50. (19) Barrio, C.; Solana, J. R. Mapping a hard-sphere fluid mixture onto a single component hard-sphere fluid. Physica A 2005, 351, 387. (20) Gidaspow, D.; Huilin, L. Equation of state and radial distribution function of FCC particles in a CFB. AIChE J. 1998, 44, 279. (21) Crowe, C.; Sommerfeld, M.; Tsuji, Y. Multiphase flows with droplets and particles; CRC Press: Boca Raton, FL, 1997. (22) Tsuji, Y.; Kawaguchi, T.; Tanaka, T. Discrete particle simulation of 2-dimensional fluidized-bed. Powder Technol. 1993, 77, 79. (23) Rhodes, M. J.; Wang, X. S.; Nguyen, M.; Stewart, P.; Liffman, K. Use of discrete element method simulation in studying fluidization characteristics, influence of interparticle force. Chem. Eng. Sci. 2001, 56, 69. (24) van Wachem, B. G. M.; van der Schaaf, J.; Schouten, J. C.; Krishna, R.; vanden Bleek, C. M. Experimental validation of Lagrangian-Eulerian simulation of fluidized beds. Powder Technol. 2001, 116, 155. (25) Feng, Y. Q.; Xu, B. H.; Zhang, S. J.; Yu, A. B.; Zulli, P. Discrete particle simulation of gas fluidization of particle mixtures. AIChE J. 2004, 50, 1713. (26) Hoomans, B. P. B.; Kuipers, J. A. M.; Briels, W. J.; van Swaaij, W. P. M. Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidized bed: a hard-sphere approach. Chem. Eng. Sci. 1996, 51, 99. (27) Ye, M.; van der Hoef, M. A.; Kuipers, J. A. M. From discrete particle model to a continuous model of Geldart A particles. Chem. Eng. Res. Des. 2005, 83, 833. (28) Li, J.; Kuipers, J. A. M. Gas-particle interactions in dense gasfluidized beds. Chem. Eng. Sci. 2003, 58, 711.

(29) Xu, B. H.; Yu, A. B. Numerical simulation of the gas-solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics. Chem. Eng. Sci. 1997, 52, 2785. (30) Helland, E.; Occelli, R.; Tadrist, L. On the heterogeneous behavior in a 2D fluidized bed. Comput. Fluid Mech. 2000, 328, 359. (31) Di Felice, R. The voidage function for fluid-particle interaction systems. Int. J. Multiphase Flow 1994, 20, 153. (32) Lun, C. K. K.; Liu, H. S. Numerical simulation of dilute turbulent gas-solid flows in horizontal channels. Int. J. Multiphase Flow 1997, 23, 575. (33) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Oxford University Press: Oxford, UK, 1980. (34) 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. (35) Jalali, P.; Polashenski, W. J.; Tynjala, T.; Zamankhan, O. Particle interactions in a dense monosized granular flow. Physica D 2002, 162, 188. (36) Huilin, L.; Shuyan, W.; Yunhua, Z.; Yang, L.; Gidaspow, D.; Ding, J. Prediction of particle motion in a two-dimensional bubbling fluidized bed using discrete hard-sphere model. Chem. Eng. Sci. 2005, 60, 3217. (37) Wen, C. Y.; Yu, Y. H. Mechanics of fluidization. Chem. Eng. Proc. Symp. Ser. 1966, 62, 100. (38) Balucani, U.; Zoppi, M. Dynamics of liquid state; Oxford Science Publications: Oxford, UK, 1994. (39) Haile, J. M.; Mansoori, J. W. Molecular-based study of fluids; American Chemical Society: Washington, DC, 1983. (40) Carnahan, N. F.; Starling, K. E. Equation of state for nonattracting rigid spheres. J. Chem. Phys. 1969, 51, 635. (41) Savage, S. B. Streaming motions in a bed of vibrationlly fluidized dry granular material. J. Fluid Mech. 1988, 194, 457. (42) Bagnold, R. A. Experiments on a gravity free dispersion of large solid spheres in a Newtonian fluid under shear. Proc. R. Soc. London 1954, A225, 49. (43) Ma, D.; Ahmadi, G. An equation of state for dense rigid sphere gases. J. Chem. Phys. 1986, 84, 3449. (44) Ghotbi, C.; Vera, J. H. Performance of three mixing rules using different equations of state for hard-spheres. Fluid Phase Equilib. 2001, 187, 321. (45) Boublik, T. Hard-sphere equation of state. J. Chem. Phys. 1970, 53, 471. (46) Mansoori, G. A.; Carnahan, K. E.; Starling, T. W.; Leland, J. Equilibrium thermodynamic properties of the mixture of hard spheres. J. Chem. Phys. 1971, 54, 1523. (47) Huilin, L.; Yurong, H.; Gidaspow, D. Hydrodynamic modeling of binary mixture in a gas bubbling fluidized bed using the kinetic theory of granular flow. Chem. Eng. Sci. 2003, 58, 1197. (48) Mathiesen, V.; Solberg, T.; Hjertager, B. H. Predictions of gasparticle flow with an Eulerian model including a realistic particle size distribution. Powder Technol. 2000, 112, 34. (49) Fan, R.; Fox, R. O. Segregation in polydisperse fluidized beds: Validation of a multi-fluid model. Chem. Eng. Sci. 2008, 63, 272.

ReceiVed for reView April 30, 2008 ReVised manuscript receiVed October 10, 2008 Accepted October 16, 2008 IE8007049