Simulation of Gas–Solid Turbulent Flow in a Cylindrical Spo

ARTICLE pubs.acs.org/EF

Computational Fluid Dynamics (CFD)
Discrete Element Method (DEM) Simulation of Gas
Solid Turbulent Flow in a Cylindrical Spouted Bed with a Conical Base Bing Ren, Wenqi Zhong,* Baosheng Jin, Zhulin Yuan, and Yong Lu School of Energy and Environment, Southeast University, Nanjing 210096, People’s Republic of China ABSTRACT: Three dimensionally coupled computational fluid dynamics (CFD) and discrete element method (DEM) were studied for modeling the turbulent gas
solid flow in a cylindrical spouted bed with a conical base. The particle motion was modeled by the DEM, and the gas motion was modeled by the k
ε two-equation turbulent model. Drag force, contact force, Saffman lift force, Magnus lift force, and gravitational force acting on an individual particle were considered in establishing the mathematical models. Calculations on the cylindrical spouted bed with an inside diameter of 152 mm, a height of 700 mm, and a conical base of 60° were carried out. Experimental results from the University of British Columbia [He, Y. L.; Qin, S. Z.; Lim, C. J.; Grace, J. R. Particle velocity profiles and solid flow patterns in spouted beds. Can. J. Chem. Eng. 1994, 72 (8), 561
568] were used as a numerical benchmark to quantitatively assess the simulations. Despite the somewhat larger simulated spout diameter found, the present simulated results were in well-agreement with the experiments. The average error of particle velocity was less than 15%. On the basis of the simulations, the development of spout with time and distributions of particle velocity, particle concentration, and spout diameters at various spouting gas velocities were obtained. Besides, detailed information on particle collision and drag forces adding on particles at different bed regions was discussed. The results showed that particle velocity gradually decreases along the radial direction, with speeds of particles moving upward decreasing with an increasing bed height, while in the annulus region, particles decelerate downward in the cylindrical section and then accelerate in the conical base. The particle concentration increases in the spout region, is kept nearly constant in the annulus region, but decreases in the fountain region along the radial axis. An increasing bed height leads to an increasing particle concentration in the spout region but a decreasing particle concentration in the annulus and fountain regions. Particle collision number, particle turbulent intensity, and the transient collision force and drag force are significantly larger in the spout region than in the annulus and fountain regions. Besides, an increasing spouting gas velocity leads to a remarkable increased speed of particles moving upward or downward, spout diameter, particle collision, drag force, and particle turbulent intensity but decreased particle concentrations.

1. INTRODUCTION Many valuable experimental and theoretical investigations on spouted beds have been carried out in the past 50 years,1 for the applications of the spouted bed for drying,2,3 granulation and coating,4,5 combustion and gasification,6
8 pyrolysis,9,10 electrochemical reaction,11,12 catalytic reaction,13,14 etc. However, the spouted bed technique has difficulty in being applied in a larger scale industrial process until now, because of some limitations, in particular the lack of full knowledge on the hydrodynamic characteristics. Numerical approaches, e.g., computational fluid dynamics (CFD) modeling, provide potentially powerful tools for realizing the complex hydrodynamics of spouted beds and the influence of the spouted bed scale, geometry, and operating conditions, as well as experiments.1 It is known that there are two main numerical approaches, i.e., the Eulerian
Eulerian approach [two-phase model (TFM)] and the Eulerian
Lagrangian approach [discrete element method (DEM)], that have been developed for understanding dense gas
solid flows. As reviewed by Zhu et al.15,16 about major applications and findings on discrete particle simulation of particulate systems, the DEM has been increasing developed and used in the past 2 decades, because r 2011 American Chemical Society

the DEM could provide detailed dynamic information at particulate dimension, e.g., the trajectories of and transient forces acting on individual particles. However, there has been little work on the DEM simulation on spouted beds; in the literature, most efforts focused on the Eulerian
Eulerian approach.17
24 Simulation of a spouted bed by the DEM is more challenging than other dense gas
solid flow systems. As Takeuchi et al.25 indicated, it is extremely difficult to establish stable spouting (i.e., stable particle circulation and entrainment26) and select a proper turbulent model of the fluid phase to describe the central turbulent spout jet. In early efforts, most are focused on two-dimensional (2D) (e.g., refs 27
30) or quasi-three-dimensional (3D) (e.g., ref 27) spouted beds. These simulations brought the unstable spouting problem, which has been confirmed by Zhao et al.28,29 through experimental and numerical simulation methods. Three-dimensional simulations on a flat spouted bed25,31 and conical-base spouted bed32 are found to successfully avoid the problem. For the model of the fluid phase, neglect of the effect of fluid turbulence would Received: June 1, 2011 Revised: July 26, 2011 Published: July 28, 2011 4095

dx.doi.org/10.1021/ef200808v | Energy Fuels 2011, 25, 4095–4105

Energy & Fuels

ARTICLE

lead to difficulties in modeling the spout
annulus interface,28 while the laminar flow model is most employed in studies, except the low Reynolds number k
ε turbulent model by Zhao et al.28,29 for a 2D spouted bed with and without a draft plate and the k
ε turbulent model by Zhong et al.26for a rectangular spout-fluid bed with a V-shaped base and Zhang et al.33 for a rectangular flat bottom spout-fluid bed. A cylindrical vessel with a conical base is the most widely used spouted bed.1 To our knowledge, there has been no published work of 3D simulation by CFD
DEM, which considered the fluid turbulence, until now. Thus, 3D CFD
DEM simulation of the gas
solid turbulent flow in a cylindrical spouted bed with a conical base is very interesting and meritorious work. The objective of this paper is to establish a 3D CFD
DEM model, which coupled a two-equation turbulent model for gas motion and a DEM model for particle motion, to simulate gas
solid turbulent flow in a cylindrical spouted bed with a conical base. Geometries of the simulated spouted bed (i.e., an inside diameter of 152 mm, a height of 700 mm, and a conical base of 60°) and the operating conditions are referenced from the experiments by the University of British Columbia.34,35 Besides, the verification of the models, development of spout with time, distributions of particle velocity, particle concentration and spout diameters at various spouting gas velocities, and detailed information on the particle collision and drag force acting on particles at different bed regions were discussed.

2. COMPUTATIONAL MODELS

The governing transport equations for turbulent energy k and turbulent dissipation rate σ are "   # ∂ μt ðεFg kÞ þ ∇ðεFg k uBg Þ ¼ ∇ ε μ þ ∇k ∂t σk þ εGk
εFg σ k þ Skd "   # ∂ μt ∇σ ðεFg σÞ þ ∇ðεFg σ uBg Þ ¼ ∇ ε μ þ σσ ∂t σ þ ε ðC1 Gk
C2 Fg σÞ þ Sσd k

Skd ¼ βjuBg
v̅ p j2 þ βðΔvp Δvp
Δug Δvp Þ Sσd ¼ C3

¼ ∂ ðεF u Þ þ ∇ðεFg uBg uBg Þ ¼
∇p þ ∇ðετÞ
Sp þ εFg gB ∂t g Bg

ð2Þ

where ε is the local void fraction, Fg is the gas density, Bug is the gas velocity vector, p is the static pressure, and Cτ is the fluid viscous stress tensor. It can be modeled as
 2 ð3Þ τ ¼ ðμ þ μt Þ ∇ uBg þ ∇ uBTg
∇ uBg I 3 where μ is the gas dynamic viscosity and μt is the turbulent viscosity μt ¼ Cμ Fg k2 =σ

ð4Þ

where Cμ is an empirically assigned constant, which is set as 0.09 in the simulation. Sp is the momentum sink given by np

Sp ¼

∑ ð BF D, i þ i¼1

ð9Þ

ðε e 0:8Þ ðε > 0:8Þ

8 > < 24 ð1 þ 0:15Rep 0:687 Þ ðRep e 1000Þ CD ¼ Rep > : 0:43 ðRep > 1000Þ

Rep ¼

εFg j uBg
B vp jdp μ

ð10Þ

ð11Þ

ð12Þ

The term β(ΔvpΔvp
ΔugΔvp) is the redistribution term, which represents the exchange of kinetic energy velocity between the particle and gas, where Δug is the gas fluctuating velocity and Δvp is the particle fluctuating velocity. According to Xiong et al.,36 the redistribution term can be rescaled as   τl τl βðΔvp Δvp
Δug Δvp Þ ¼
2βk 1
ð13Þ τ l þ τ d τl þ τd where τl is the Lagrangian time scale of the gas phase and τd is the response time scale of the particle phase. They are calculated by τd ¼

ð5Þ

where np is the number of particles in a considered computational cell FLM,i, and B FLS,i are and ΔV is the volume of the computational cell. B FD,i, B the drag force, Magnus lift force, and Saffiman lift force acting on particle i, respectively.

σσ k S k d

8 μð1
εÞ > > > < εd 2 ½150ð1
εÞ þ 1:75Rep  p β¼ > 3 μð1
εÞ
2:7 > > : 4 CD d 2 ε Rep p

B F LS, i þ B F LM, i Þ ΔV

ð8Þ

2 The term β|u Bg
Bvp| is the generation term caused by particle resistance. Bvp is the mean velocity of the particles in the unit volume. β is the gas
solid interphase drag coefficient, defined as

motion are used to compute the motion of the gas phase, given respectively by ð1Þ

ð7Þ

In these equations, Gk represents the generation of turbulence kinetic energy because of the mean velocity gradients, C1 = 1.44 and C2 = 1.92, and σk = 1.0 and σσ = 1.3 are the turbulent Prandtl numbers for k and σ, respectively. Skd and Sσd represent the influence of the particle to the fluid, calculated as described in eqs 8 and 9.

2.1. Gas Phase. The local averaged equations of continuity and

∂ ðεFg Þ þ ∇ðεFg uBg Þ ¼ 0 ∂t

ð6Þ

4dp Fp 3CD Fg j uBg
B v pj

τl ¼ 0:35

k σ

ð14Þ

ð15Þ

C3 in eq 9 is an empirically assigned constant, and here, C3 = 1.2, cited from the channel flow experiment by Pourahmadi and Humphery.37 4096

dx.doi.org/10.1021/ef200808v |Energy Fuels 2011, 25, 4095–4105

Energy & Fuels

ARTICLE

2.2. Particle Phase. According to Newton’s equation for motion, the translational and rotational motions of particle i are given by mi

Ii

dB v pi dt

¼ mi gB þ B F C, i þ B F D, i þ B F LS, i þ B F LM, i

dw Bpi ¼ dt

ni

∑ ð TBtij þ j¼1

T Bnij Þ

ð16Þ

B F cnij

Req

1 1 ¼ þ Ri Rj

Sn ¼ 2Eeq

B F ctij ¼ St δtij

B F C, i ¼

ð28Þ

n

∑ ð BF cnij þ j¼1

B F ctij þ B F dnij þ B F dtij Þ

ð20Þ

B F D, i ¼

β ðu
v Þ Fg Bg Bp, i

ð30Þ

where Bug, Bvp,i , and Fg are fluid velocity, velocity of particle i, and fluid density, respectively. 2.2.3. Saffman Lift Force. The Saffman lift force is calculated by ffi vffiffiffiffiffiffiffiffiffi u    u∂ug  ∂ug 0:5 2 vp, i ÞðFg μg Þ dp CLS t sgn B F LS, i ¼ 1:615ðuBg
B  ∂n  ∂n ð31Þ in which the Saffman lift coefficient CLS is calculated from the correlation proposed by Mei42 8 < ð1
0:3314ς0:5 Þexpð
0:1Rep Þ þ 0:3314ς0:5 ðRep e 40Þ CLS ¼ : 0:0524ς0:5 ðRep Þ0:5 ðRep > 40Þ ð32Þ

ð21Þ where

qffiffiffiffiffiffiffiffiffiffiffiffiffi St ¼ 8Geq Req δnij

Geq is the equivalent shear modulus, defined as !
1 1
γj 1
γi þ Geq ¼ Gi Gj

The tangential microslip is neglected. The gross sliding condition provided by Coulomb’s law of friction is imposed as a constraint to Bcnij|, the tangential force is the tangential force calculation. If |F Bctij| > μ|F given by v tij B jB v tij j

ð33Þ

w v 1 B F LM, i ¼ Fg vr 2 πdp 2 CLM Br Br jw v rj 8 Br jj B

ð34Þ

where vr and wr are the relative velocity and relative rotation angular velocity between the gas and particle, respectively. CLM is the Magnus lift, which is calculated by

CLM

ð26Þ

  ∂u 0:5dp   ¼   jð uBg
B v p, i Þj∂n

2.2.4. Magnus Lift Force. The Magnus lift force is calculated by the following correlation:43,44

ð24Þ

and the correspondent value of the tangential stiffness tends to 0. The normal and tangential damping forces are sited from Raji,41 which can be described as rffiffiffi 5 pffiffiffiffiffiffiffiffiffiffiffi B F dnij ¼
2 ð25Þ v nij ψ Sn meq B 6 rffiffiffi 5 pffiffiffiffiffiffiffiffiffiffiffi ψ St meq B v tij 6

CLS

ð22Þ

ð23Þ

ð29Þ

2.2.2. Drag Force. The drag force for a single particle is calculated as

with

B F dtij ¼
2

qffiffiffiffiffiffiffiffiffiffiffiffiffi Req δnij

#
1

with Ei, γi, and Ri and Ej, γj, and Rj being the Young’s modulus, Poisson ratio, and radius of each particle in contact. The tangential contact force B Fctij depends upon the tangential overlap δtij and the tangential stiffness St

F cnij j B F ctij ¼ μj B

ð27Þ

where e is the coefficient of restitution. The total force acting on particle i because of particle collision is calculated as a combination of the contact force and damping force, which is expressed as

ð18Þ

where the equivalent Young’s modulus Eeq and equivalent radius Req are defined as " #
1 ð1
γj 2 Þ ð1
γi 2 Þ þ ð19Þ Eeq ¼ Ei Ej "

ln e ψ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ln e þ π2

ð17Þ

where B FC,i is the contact force, B FD,i is the drag force, B FLS,i is the Saffman lift force, and B FLM,i is the Magnus lift force.Bvpi is the velocity of particle i. wpi is the rotational velocity of particle i. Ii is the inertia of the particle. B Tnij are generated by the tangential forces and the Torques B Ttij and B rolling friction, respectively.38 2.2.1. Contact Force. The Hertz
Mindlin no-slip model is employed in our current study. It is based on the classical Hertz’s theory39 for the normal direction and simplifications of the model developed by Mindlin and Deresiewicz40 for the tangential direction. The normal contact force, B Fcnij, is a function of the normal overlap δnij and is given by 4 pffiffiffiffiffiffiffi ¼ Eeq Req δnij 3=2 3

where meq = ((1/mi) + (1/mj))
1 is the equivalent mass and vnij and vtij are the relative normal and tangential velocities, respectively. ψ and Sn are the damping ratio and the normal stiffness, given by

8 > Br j >j w > < j v j dp B r ¼ jw > B rj
0:522 > > Þ : j v j dp ð0:178 þ 0:822Rep Br

ðRep e 1Þ ð1 < Rep < 1000Þ

ð35Þ

3. COMPUTATIONAL CONDITIONS The simulated vessel has a similar geometry as the previous experiments taken by He et al.;34,35 it is a fully cylindrical Plexiglas column with an inside diameter of 152 mm and a height of 700 mm with a 60° conical base, shown in Figure 1a. 4097

dx.doi.org/10.1021/ef200808v |Energy Fuels 2011, 25, 4095–4105

Energy & Fuels

ARTICLE

Figure 1. Geometry of the vessel and numerical grids: (a) sketch of the simulated bed, (b) numerical grid for the gas, and (c) numerical grid for the particle.

Table 1. Physical and Numerical Parameters properties

value

particle diameter, dp (mm)

2.9

particle density, Fp (kg/m3)

1450

gas density, Fg (kg/m3)

1.205

gas viscosity, μg (Pa/s)

18.1  10
6

initial bed height, H0 (mm)

325

vessel height, H (mm)

700

diameter of the bed, Dt (mm) diameter of the spout gas inlet, Ds (mm)

152 19

minimum spouting velocity, ums (m/s)

1.15

spouting gas velocity, us (m/s)

1.1
1.3ums

number of particles, N

217552

restitution conefficient, e

0.9

are checked on the basis of the cells that contain two or more particles. If there is no more than one particle in each grid cell, then no contact detection needs to take place. By this way, the speed of the simulation is significantly improved. In this work, the length of the grid cells is set to 4dp. Detailed physical and numerical parameters chosen in this simulation are listed in Table 1. The number of the particles corresponding to the same initial bed height of 325 mm as in the experiments is about 2 000 000. To reduce the time of simulation, we reduce the number of particles in our simulation by assuming several closely spaced particles as a composite spherical particle. The radius and density of the new particle can be calculated as follows:

friction coefficient particle
particle, μp paritcle
wall, μw Possion’s ratio

0.3 0.25

particle, γp

0.3

wall, γw

0.3

shear modulus particle, Gp (N/m2) 2

wall, Gw (N/m ) time step of calculation, Δt (s)

3  109 3  109 1  106

Figure 1b shows the gird used in the fluid field. Three grid domains were tested, containing 25 568, 69 510, and 99 960 CFD cells, respectively. The difference is less than 3% for the simulation results; thus, the computational domain containing 69 510 CFD cells is chosen in this work. The main computational challenge in DEM simulation is the detection of contacts. In the DEM procedure, the domain is divided into uniform grid cells, as shown in Figure 1c. Contacts

1 1 Na πdp0 3 ¼ πdp 3 ð1
εÞ 6 6

ð36Þ

1 1 Na πdp0 3 Fp0 ¼ πdp 3 Fp 6 6

ð37Þ

where Na is the hypothetical number of assembling particles, and in this simulation, the number is set to 5, to make this hypothesis both cost-effective and practical. dp0 and Fp0 are the initial particle diameter and density, respectively, and in the experiment by He et al.,34,35 dp0 = 1.4 mm and Fp0 = 2503 kg/m3. The voidage ε is 0.42 in our packing test. Then, we obtain the radius and density of the new particle, dp = 2.9 mm and Fp = 1450 kg/m3. The minimum spouting velocity based on our simulation is 1.15 m/s; it is calculated through the equation proposed by Choi and Meisen.45 !0:263 Fp
Fg ums ¼ 18:5 ðH0 =Dt Þ
0:103 ðdp =Dt Þ1:19 Fg pffiffiffiffiffiffiffiffiffiffi ð38Þ ðDs =Dt Þ0:373 2gH0 Applicability of the correlation has been examined by predicting the minimum spouting velocity of glass beads with a diameter of 4098

dx.doi.org/10.1021/ef200808v |Energy Fuels 2011, 25, 4095–4105

Energy & Fuels

ARTICLE

1.41 mm and a density of 2503 kg/m3 in the experiment by He et al.;34,35 the value obtained in the correlation is 0.56 m/s for the static bed height of 325 mm, which agrees well with the experimental value of 0.54 m/s. The calculations were made for gas velocities of 1.1ums, 1.2ums, and 1.3ums, to compare the results to the experiment. The time step Δt is determined from the viewpoint of the Rayleigh wave speed of force transmission on the surface of elastic bodies. It is assumed that, during a single time step, the disturbance can not propagate from the particle farther than its immediate neighbor particles. In accordance with the assumption, the time step should be smaller than the critical Rayleigh time step TR, which is given by  1=2 Fp 1 TR ¼ πdp ð0:1631γ þ 0:8766Þ
1 2 G

ð39Þ

In the simulation, the time step is set to 1  10
6 s, about 20% of the critical Rayleigh time step. The calculation of simulations starts with a packed bed with the voidage of about 0.42 and the height of 325 mm. The fully developed gas is introduced into the bottom of the packed bed. On the inlet side, a solid plane wall normal to the z axis was assumed, except for the opening of the nozzle. The no-slip condition was applied at the wall. For the outflow boundary, Sommerfeld’s non-reflective condition and Neuman’s condition ∂p/∂z = 0 were applied.26 In this work, the gas phase was solved by the non-staggered semi-implicit method for pressure-linked equations (SIMPLE) method46 and the time advancement of pressure and velocity field is based on the simplified marker and cell (SMAC) method,47 which has been successfully applied to the simulation of multiphase turbulence flow. The translational and rotational motions of discrete particles are solved by the explicit time integration. The coupling scheme between DEM and CFD can be described as follows: At the beginning, the flow field of the gas phase is resolved by the CFD solver. When the transient calculation is iterated to convergence for a time step, the flow field is passed to the coupling module. The particle
fluid forces acting on each particle are calculated, and the particle positions will be updated. The new particle positions are passed back to the coupling module, and the momentum exchange is added to each of the mesh cells to represent the effect of the interphase interaction between the solid and the gas per unit volume.

4. RESULTS AND DISCUSSION 4.1. Flow Patterns in the Spouted Bed. A full development of a spout in the spouted bed at us = 1.3ums is presented in Figure 2. It is based on a total of a 3 s simulation, and a series of instantaneous snapshots are generated every 0.001 s. Different colors are used to distinguish the different magnitude of particle velocity. The color red represents the maximum particle velocity, while the color blue represents the minimum particle velocity. The pictures are obtained by projecting particles on a 3 mm depth cross-section in the center of the spouted bed. Figure 2a is the initial packed bed; fully developed gas is introduced into the bottom of the packed bed at this moment. Particles are transported upward by the spouting gas; an obvious jet will occur inside the bed, as seen in Figure 2b. The internal spout expands upward as in Figure 2c and then breaks through the peripheral bed level. A fountain zone is formed in Figure 2d; it grows

Figure 2. Development of a spout in the spout bed (us = 1.3ums): (a) t = 0 s, (b) t = 0.1 s, (c) t = 0.2 s, (d) t = 0.3 s, (e) t = 0.4 s, (f) t = 0.5 s, (g) t = 0.6 s, (h) t = 0.7 s, (i) t = 0.8 s, (j) t = 0.9 s, (k) t = 1 s, (l) t = 1.1 s, (m) t = 1.2 s, (n) t = 1.3 s, (o) t = 1.5 s, and (p) t = 2.5 s.

gradually with time, as present from panels e to n of Figure 2. Finally, the spout reaches a stable and nonpulsating state, as in panels o and p of Figure 2. The flow area can be obviously divided into three regions: a central spout region, where the gas and particles rise at high velocity and the particle concentration is low, a fountain zone, where particles rise to their highest positions and then rain back onto the surface of the annulus, and an annulus zone between the spout and the column wall, where particles move slowly downward as a dense phase, with counter-current percolation of the fluid. 4.2. Radial Distribution of Particle Velocity. Figure 3 presents the time-averaged vector field of the particle velocity of a fully developed spout at us = 1.3ums. A systematic cyclic of solid movements exists in the spout, annulus, and fountain regions. It can be seen that particles from the annulus feed into the spout over the entire spout length, especially near the conical base, which agrees well with the previous experimental investigation in the conventional spouted bed34,48,49 and also in the conical spouted bed.50,51 While on the top of the fountain region, particles are cascading downward. Detailed analyses of the particle velocity are discussed below with the help of Figures 4
6. Comparisons of simulated radial profiles of time-averaged vertical particle velocities to experiments at us = 1.3ums are shown in panels a
c of Figure 4. The simulation results agree well with the experiments. In the spout region, the error is larger at the center and increases with the bed height. In the annular, the predicted particle velocities are somewhat larger than the experimental data, which might be caused by the increase of the particle 4099

dx.doi.org/10.1021/ef200808v |Energy Fuels 2011, 25, 4095–4105

Energy & Fuels

ARTICLE

Figure 3. Time-averaged vector field of the particle velocity of a fully developed spout (us = 1.3ums and t = 3 s).

diameter. In the fountain region, the simulated particle velocities almost coincide entirely with the experimental data, except in the central region at the height of 370 mm. The mean relative errors of predictions to the experimental data are 10.9, 9.8, and 13.4% in the spout, annulus, and fountain regions, respectively. The overall mean relative error of simulated particle velocities is 12.1%. The comparisons show that the present models are helpful for computing the right trends of the complex gas
solid flow in the spouted bed. It can also be clearly seen in panels a
c of Figure 4 the axial and radial distributions of the particle velocity in the three regions of the spouted bed. It is shown that the particle velocity decreases along the radial direction in the whole bed. The decreasing trend is more remarkable in the spout region than in the fountain region, while in the annulus region, the speed of particles moving downward is much smaller. In the axial direction, the distribution of the particle velocity is different in the three regions. In the spout region, with the increase of the bed height, particle velocity decreases along the central axis but increases at the periphery of the spout. In the annular region, particles decelerate downward in the cylindrical section, while they accelerate when traveling into the conical base. While in the fountain region, both speeds of particle move upward and downward, decreasing with bed elevation. The trend of the radial velocity profile in the fountain region is similar to the quasi-3D spouted bed by Kawaguchi et al.27 and the flat bottom spouted bed by Takeuchi et al.25 Figures 5 and 6 show the comparison of the particle velocity profiles at various gas velocities at bed heights of 268 and 400 mm, respectively.An increasing spouting gas velocity leads to a

Figure 4. Comparison of simulated particle velocities to experimental results in the (a) spout region, (b) annulus region, and (c) fountain region.

remarkable increased speed of particles moving upward and downward throughout the entire spouted bed. The overall particle velocity profiles in the spout and annulus regions at the height of 268 mm for different gas velocities are quite similar in shape in Figure 5. Particle velocity increases from 2.38 to 3.35 m/s with gas velocity increasing from 1.1ums to 1.3ums in the center of the bed. The changes of particle velocity in the annulus are much smaller, however, by zooming in on the curves. It is obviously seen that the speed of the particle downward increases from 0.01 to 0.02 m/s when the gas velocity increases from 1.1ums to 1.3ums by zooming in on the curves. The changes of particle velocity with gas velocity in the fountain region is presented in Figure 6, with spouting gas velocity increasing from 4100

dx.doi.org/10.1021/ef200808v |Energy Fuels 2011, 25, 4095–4105

Energy & Fuels

ARTICLE

Figure 5. Simulation result of the radial distribution of the particle velocity at various gas velocities (z = 268 mm).

Figure 8. Radial distribution of the particle concentration at various gas velocities (z = 268 mm).

Figure 6. Simulation result of the radial distribution of the particle velocity at various gas velocities (z = 400 mm).

Figure 9. Radial distribution of the particle concentration at various gas velocities (z = 400 mm).

Figure 7. Radial distributions of the particle concentration at different bed heights (us = 1.3ums).

1.1ums to 1.3ums, particle velocity increasing from 0.74 to 1.68 m/s in the fountain core, and in the fountain outer downflow region, the speed of the particle increasing from 0.31 to 1.12 m/s. 4.3. Radial Distribution of the Particle Concentration. The simulation values of the particle concentration at us = 1.3ums at different bed heights are shown in Figure 7. The particle concentration in the spout region at z = 318 and 118 mm are

axisymmetric, similar to the parabolic shapes, and it increases with an increasing radial distance. In the interface of the spout and annulus, the particle concentration reaches the peak value and then slightly decreases. In the annulus region, the particle concentration fluctuates in a very small range. At z = 420 and 620 mm in the fountain region, the particle concentration is higher in the core region and decreases gradually with an increasing radial distance. In the axial direction, the particle concentration gradually increases with height in the spout region but decreases in the annulus and fountain regions. In previous studies, some authors have observed that the maximum particle concentration is in the upper surface near the wall. However, the location of the maximum particle concentration is not very clearly seen in the present simulations, because it is sensitive to many factors, i.e., operating condition, particle diameter, and geometry of the bed. Comparisons of radial profiles of the particle concentration for three different gas velocities, us/ums = 1.1, 1.2, and 1.3, are presented in Figures 8 and 9. An increasing gas velocity leads to the decrease of the particle concentration throughout the entire bed. The particle concentration increases uniformity with a decreasing gas velocity over the spout and annulus regions, while in the fountain region, the particle concentration increases more obviously in the core than in the outer downflow region. It indicated that the increasing gas velocity improves the gas
solid movement over the entire section. 4101

dx.doi.org/10.1021/ef200808v |Energy Fuels 2011, 25, 4095–4105

Energy & Fuels

Figure 10. Comparisons of simulated spout diameter to experimental results.

ARTICLE

Figure 12. Radial distribution of the collision number (z = 350
400 mm).

Figure 11. Radial distribution of the collision number (z = 150
200 mm).

4.4. Spout Diameter. The spout diameter can be considered as one of the key parameters to split the flow distribution between the spout and annulus regions. In the present simulation, the spout diameter can be determined by tracing the particle velocity over the entire spout and annulus regions and detecting the point where particle velocities go through zero. The corresponding radial positions at different bed levels therefore represent the boundary between the spout and annulus regions. The spout borders where the vertical particle velocity is zero are outlined in the lower left of Figure 10, with the help of a steady spout condition snapshot at us/ums = 1.3 and t = 3 s. A comparison of experimental results by He et al.35 to the predicted results by Kawaguchi et al.27 and the present simulation is shown in Figure 10. A qualitative difference exists between both the calculation and the experiments. The spout boundary for each of the three gas flow rates in the experiment diverges in the conical base region, narrows further up the column, and then diverges again near the bed surface. The simulated spout shape by Kawaguchi et al.27 is different from the experiments, which they call “reverse conical shaped”. They explained that this discrepancy was due to the difference of the particle size. The present discrepancy also seemed to be caused by this. The spout diameter increases monotonously with the height. The increment is less than that from the simulation by Kawaguchi et al.27 but is still different from the nearly constant value at the bed

Figure 13. Average collision force acting on a single particle: (a) z = 150
200 mm and (b) z = 350
400 mm.

height between 100 and 250 mm in the experiment. However, the trend of the spout diameter for different gas velocities agrees well with the experiments in both simulations; that is, the spout diameter increases with the increasing gas velocity. 4.5. Collision Number. The collision number is the number of completed collisions between particles in one time step. Radial distributions of the collision number at different gas velocities are plotted in Figures 11 and 12; the sample section is at the height of 150
200 and 350
400 mm, respectively. An increasing gas velocity significantly increases the collision number at the spout 4102

dx.doi.org/10.1021/ef200808v |Energy Fuels 2011, 25, 4095–4105

Energy & Fuels

Figure 14. Average drag force acting on a single particle: (a) z = 150
200 mm and (b) z = 350
400 mm.

and annulus regions (z = 150
200 mm) and the fountain region (z = 350
400 mm). The radial distributions of the collision number have different characteristics in different bed elevations. At the height of 150
200 mm, the collision number is much larger in the center of the spouted bed. This indicated that particles collide with each other much more intensely in the spout region than in the annulus region. The peak of the collision number exists near the boundary between the spout and annulus regions, probably because of the strong entrainment of particles from the annulus to the spout region. However, the collision number does not vary much at the height of 350
400 mm. At us = 1.3ums, it is slightly higher in the bed center than near the wall. While at us = 1.2ums and 1.1ums, the collision number changes are small in the whole fountain region; the average value is 7.6 and 3.4, respectively. The result indicates that, in the fountain region, the collision between particles becomes weakened and tends to be uniform along the radial axis. 4.6. Forces Acting on Particles. Figure 13 presents the dimensionless ratio of the normal collision force Fcn to the gravity of the individual particle Gp at different bed heights. It can be seen that the fluctuation of the ratio is much larger in the central region of the bed at the height of 150
200 mm, which indicted that the collision in the spout region is more intense and disorderly than that in the annulus region. While at the height of 350
400 mm, the ratio fluctuates along its average value of
0.011 and the variation is about 20 orders of magnitude lower than that in the spout region. This reveals that the collision between particles tends to be mild in the fountain region and is homogeneous along the radial direction.

ARTICLE

Figure 15. Radial distribution of the particle turbulent intensity at various vessel heights: (a) z = 200 mm and (b) z = 400 mm.

As shown in Figure 14, the dimensionless ratio of the average normal drag force FDn to gravity Gp shares the same characteristics. The fluctuation of the drag force is significantly larger in the spout region than in the annulus and fountain regions. It should be noted that the changes of the drag force are larger than that of the collision force in the spout region, which indicated that the turbulent gas
solid flow is stronger than the collision between particles. 4.7. Particle Turbulent Intensity. The particle turbulent intensity radial profiles at a bed height of 200 and 400 mm are shown in Figure 15. The particle turbulent intensity is defined as ((Δvx2 + Δvy2 + Δvz2)/3)1/2/vp, in which Δvx, Δvy, and Δvz are the changes of the particle velocity in one time step along the direction x, y, and z and vp is the mean particle velocity. For the spout and annulus regions (z = 200 mm), particle turbulent intensity is always much higher in the spout region and the interface between the spout and annulus regions, while it tends to be constant in the annulus region. For the fountain region (z = 400 mm), the particle turbulent intensity fluctuates slightly along the radial direction. For all of the three regions, an increasing gas velocity leads to a larger particle intensity. The simulation results indicate that the particle
particle turbulent interaction is larger in the spout region than that in the fountain and annulus regions and becomes more intense with an increasing gas velocity.

5. CONCLUSION A three-dimensional CFD
DEM model was established to simulate gas
solid turbulent flow in a cylindrical spouted bed with a conical base. The particle motion was modeled by the 4103

dx.doi.org/10.1021/ef200808v |Energy Fuels 2011, 25, 4095–4105

Energy & Fuels DEM, and the gas motion was modeled by the k
ε two-equation turbulent model. The drag force, contact force, Saffman lift force, Magnus lift force, and gravitational force acting on individual particles were considered in establishing the mathematical models. The present simulated results were in well-agreement with published experimental results. The following conclusions can be drawn on the basis of the simulations: (1) The particle velocity decreases gradually along the radial direction. Speeds of particle moving upward in the spout and fountain regions decrease with bed elevation. In the annular region, particles decelerate downward in the cylindrical section, while they accelerate when traveling into the conical base. (2) The particle concentration increases along the radial direction in the spout region but decreases in the fountain region, while it is nearly kept constant in the annulus region. With bed elevation, the particle concentration gradually increases in the spout region but decreases in the annulus and fountain regions. (3) Particle collisions and turbulent intensity in the spout region are more intense than those in the annulus and fountain regions. Besides, both the fluctuation transient collision force and drag force are significantly larger in the spout region than in the annulus and fountain regions. (4) An increasing spouting gas velocity leads to the remarkably increased speed of particles moving upward or downward, spout diameter, particle collision, particle turbulent intensity, transient collision force, and drag force but decreased particle concentration.

’ AUTHOR INFORMATION Corresponding Author

*Telephone: +86-25-83794744. Fax: +86-25-83795508. E-mail: [email protected].

’ ACKNOWLEDGMENT Financial support from the National Natural Science Foundation of China (50976025), the Major State Basic Research Development Program of China (2011CB201505), and the Scientific Research Foundation of Graduate School of Southeast University (YBJJ1107) is sincerely acknowledged. ’ NOMENCLATURE CLM = Magnus lift coefficient CLS = Saffman lift coefficient D = diameter of the bed (mm) Dt = diameter of the spout gas inlet (mm) dp = particle diameter (mm) Eeq = equivalent Young’s modulus (N/m2) FC,i = contact force (N) FD,i = drag force (N) FLS,i = Saffman lift force (N) FLM,i = Magnus lift force (N) Fcnij = normal contact force (N) Fctij = tangential contact force (N) Fdnij = normal damping force (N) Fdtij = tangential damping force (N) e = restitution coefficient Geq = equivalent shear modulus (N/m2) H0 = initial bed height (mm) H = vessel height (mm) Ii = particle motion of inertia (N m) k = turbulent energy (m2/s2)

ARTICLE

mi = mass of particle i (kg) meq = equivalent mass (kg) N = number of particles Na = number of assembling particles Nc = number of collisions np = particle number in a computational cell ni = number of particles in contact with particle i p = pressure drop (Pa) Req = equivalent particle contact radius (m) Rep = particle Reynolds number r = local position in the radial direction (mm) rs = radius of the spout region (mm) Sp = momentum sink (N m
3) St = tangential stiffness (N m) Sn = normal stiffness (N m) TR = critical Rayleigh time step (s) ug = gas velocity (m/s) ums = minimum spouting velocity (m/s) us = spouting gas velocity (m/s) vp = particle velocity (m/s) vnij = relative normal velocity between particle i and particle j (m/s) vtij = relative tangential velocity between particle i and particle j (m/s) z = bed height (mm) γp = particle Poisson’s ratio γw = wall Poisson’s ratio δnij = normal displacements between particle i and particle j (m) δtij = tangential displacements between particle i and particle j (m) ψ = damping ratio Δug = gas fluctuating velocity (m/s) Δvp = particle fluctuating velocity (m/s) Δt = time step of calculations (s) μ = gas dynamic viscosity (Pa s) μt = gas turbulent viscosity (Pa s) μg = gas dynamic viscosity (Pa s) μp = particle
particle friction coefficient μw = particle
wall friction coefficient ε = void fraction Fg = gas density (kg/m3) Fp = particle density (kg/m3) τl = Lagrangian time scale of the gas phase τd = response time scale of the particle phase τ = gas turbulence stress (Pa) σ = turbulent dissipation rate (m2/s3)

’ REFERENCES (1) Epstein, N.; Grace, J. R. Spouted and Spout-Fluid Beds: Fundamentals and Applications; Cambridge University Press: New York, 2011. (2) Mathur, K. B.; Gishler, P. E. A study of the application of the spouted bed technique to wheat drying. J. Appl. Chem. 1955, 5 (11), 624–636. (3) Zahed, A. H.; Epstein, N. Batch and continuous spouted bed drying of cereal grains: The thermal equilibrium model. Can. J. Chem. Eng. 1992, 70 (5), 945–953. (4) Liu, L. X.; Litster, J. D. Coating mass distribution from a spouted bed seed coater: Experimental and modelling studies. Powder Technol. 1993, 74 (3), 259–270. (5) Jono, K.; Ichikawa, H.; Miyamoto, M.; Fukumori, Y. A review of particulate design for pharmaceutical powders and their production by spouted bed coating. Powder Technol. 2000, 113 (3), 269–277. (6) Zhao, J.; Lim, C. J.; Grace, J. R. Flow regimes and combustion behaviour in coal-burning spouted and spout-fluid beds. Chem. Eng. Sci. 1987, 42 (12), 2865–2875. 4104

dx.doi.org/10.1021/ef200808v |Energy Fuels 2011, 25, 4095–4105

Energy & Fuels (7) Lim, C. J.; Watkinson, A. P.; Khoe, G. K.; Low, S.; Epstein, N.; Grace, J. R. Spouted, fluidized and spout-fluid bed combustion of bituminous coals. Fuel 1988, 67 (9), 1211–1217. (8) Abdul Salam, P.; Bhattacharya, S. C. A comparative study of charcoal gasification in two types of spouted bed reactors. Energy 2006, 31 (2
3), 228–243. (9) Lisb^oa, A. C. L.; Watkinson, A. P. Pyrolysis with partial combustion of oil shale fines in a spouted bed. Can. J. Chem. Eng. 1992, 70 (5), 983–990. (10) Aguado, R.; Olazar, M.; San Jose, M. J.; Aguirre, G.; Bilbao, J. Pyrolysis of sawdust in a conical spouted bed reactor. Yields and product composition. Ind. Eng. Chem. Res. 2000, 39 (6), 1925–1933. (11) Hadzismajlovic, D. E.; Pavlovic, M. G.; Popov, K. I. The annulus of a spouted bed as a three-dimensional electrode. Hydrometallurgy 1989, 22 (3), 393–401. (12) Jiricny, V.; Roy, A.; Evans, J. Copper electrowinning using spouted-bed electrodes: Part I. Experiments with oxygen evolution or matte oxidation at the anode. Metall. Mater. Trans. B 2002, 33 (5), 669–676. (13) Mathur, K. B.; Lim, C. J. Vapour phase chemical reaction in spouted beds: A theoretical model. Chem. Eng. Sci. 1974, 29 (3), 789–797. (14) Rovero, G.; Piccinini, N.; Grace, J. R.; Epstein, N.; Brereton, C. M. H. Gas phase solid-catalysed chemical reaction in spouted beds. Chem. Eng. Sci. 1983, 38 (4), 557–566. (15) Zhu, H. P.; Zhou, Z. Y.; Yang, R. Y.; Yu, A. B. Discrete particle simulation of particulate systems: Theoretical developments. Chem. Eng. Sci. 2007, 62 (13), 3378–3396. (16) Zhu, H. P.; Zhou, Z. Y.; Yang, R. Y.; Yu, A. B. Discrete particle simulation of particulate systems: A review of major applications and findings. Chem. Eng. Sci. 2008, 63 (23), 5728–5770. (17) Huilin, L.; Yurong, H.; Wentie, L.; Ding, J.; Gidaspow, D.; Bouillard, J. Computer simulations of gas
solid flow in spouted beds using kinetic-frictional stress model of granular flow. Chem. Eng. Sci. 2004, 59 (4), 865–878. (18) Du, W.; Bao, X.; Xu, J.; Wei, W. Computational fluid dynamics (CFD) modeling of spouted bed: Influence of frictional stress, maximum packing limit and coefficient of restitution of particles. Chem. Eng. Sci. 2006, 61 (14), 4558–4570. (19) Shirvanian, P. A.; Calo, J. M.; Hradil, G. Numerical simulation of fluid
particle hydrodynamics in a rectangular spouted vessel. Int. J. Multiphase Flow 2006, 32 (6), 739–753. (20) Wu, Z. H.; Mujumdar, A. S. CFD modeling of the gas
particle flow behavior in spouted beds. Powder Technol. 2008, 183 (2), 260–272. (21) Bettega, R.; de Almeida, A. R. F.; Corr^ea, R. G.; Freire, J. T. CFD modelling of a semi-cylindrical spouted bed: Numerical simulation and experimental verification. Can. J. Chem. Eng. 2009, 87 (2), 177–184. (22) Gryczka, O.; Heinrich, S.; Deen, N. G.; van Sint Annaland, M.; Kuipers, J. A. M.; M€orl, L. CFD modeling of a prismatic spouted bed with two adjustable gas inlets. Can. J. Chem. Eng. 2009, 87 (2), 318–328. (23) Duarte, C. R.; Olazar, M.; Murata, V. V.; Barrozo, M. A. S. Numerical simulation and experimental study of fluid
particle flows in a spouted bed. Powder Technol. 2009, 188 (3), 195–205. (24) Santos, K. G.; Murata, V. V.; Barrozo, M. A. S. Threedimensional computational fluid dynamics modelling of spouted bed. Can. J. Chem. Eng. 2009, 87 (2), 211–219. (25) Takeuchi, S.; Wang, S.; Rhodes, M. Discrete element simulation of a flat-bottomed spouted bed in the 3-D cylindrical coordinate system. Chem. Eng. Sci. 2004, 59 (17), 3495–3504. (26) Zhong, W. Q.; Xiong, Y. Q.; Yuan, Z. L.; Zhang, M. Y. DEM simulation of gas
solid flow behaviors in spout-fluid bed. Chem. Eng. Sci. 2006, 61 (5), 1571–1584. (27) Kawaguchi, T.; Sakamoto, M.; Tanaka, T.; Tsuji, Y. Quasithree-dimensional numerical simulation of spouted beds in cylinder. Powder Technol. 2000, 109 (1
3), 3–12. (28) Zhao, X. L.; Li, S. Q.; Liu, G. Q.; Yao, Q.; Marshall, J. S. DEM simulation of the particle dynamics in two-dimensional spouted beds. Powder Technol. 2008, 184 (2), 205–213.

ARTICLE

(29) Zhao, X. L.; Li, S. Q.; Liu, G. Q.; Song, Q.; Yao, Q. Flow patterns of solids in a two-dimensional spouted bed with draft plates: PIV measurement and DEM simulations. Powder Technol. 2008, 183 (1), 79–87. (30) Rong, L.; Zhan, J. Improved DEM
CFD model and validation: A conical-base spouted bed simulation study. J. Hydrodyn., Ser. B 2010, 22 (3), 351–359. (31) Takeuchi, S.; Shan Wang, X.; Rhodes, M. J. Discrete element study of particle circulation in a 3-D spouted bed. Chem. Eng. Sci. 2005, 60 (5), 1267–1276. (32) Takeuchi, S.; Wang, S.; Rhodes, M. Discrete element method simulation of three-dimensional conical-base spouted beds. Powder Technol. 2008, 184 (2), 141–150. (33) Zhang, Y.; Jin, B.; Zhong, W.; Ren, B.; Xiao, R. DEM simulation of particle mixing in flat-bottom spout-fluid bed. Chem. Eng. Res. Des. 2010, 88 (5
6), 757–771. (34) He, Y. L.; Lim, C. J.; Grace, J. R.; Zhu, J. X.; Qzn, S. Z. Measurements of voidage profiles in spouted beds. Can. J. Chem. Eng. 1994, 72 (2), 229–234. (35) He, Y. L.; Qin, S. Z.; Lim, C. J.; Grace, J. R. Particle velocity profiles and solid flow patterns in spouted beds. Can. J. Chem. Eng. 1994, 72 (4), 561–568. (36) Xiong, Y. Q.; Yuan, Z. L.; Zhang, M. Y. Three dimensional numerical simulation on conveying properties of gas
solid injector under pressurization. Chin. J. Chem. Ind. Eng. 2004, 55, 1638–1642. (37) Pourahmadi, F.; Humphery, J. A. C. Prediction of curved channel flow with k
ε model of turbulence. AIAA J. 1983, 7, 59–75. (38) Kuang, S. B.; Yu, A. B.; Zou, Z. S. Computational study of flow regimes in vertical pneumatic conveying. Ind. Eng. Chem. Res. 2009 48 (14), 6846–6858. (39) Hertz, H. Uber die beruhrung fester elastischer korper (On the contact of elastic solids). J. Reine Angew. Math. 1882, 92, 156–171. (40) Mindlin, R. D.; Deresiewicz, H. Elastic spheres in contact under varying oblique forces. J. Appl. Mech. 1953, 20, 327. (41) Raji, A. O. Discrete element modelling of the deformation of bulk agricultural particulates. Ph.D. Thesis, Department of Agricultural and Environmental Sciences, University of Newcastle, Newcastle upon Tyne, U.K., 1999. (42) Mei, R. An approximate expression for the shear lift force on a spherical particle at finite Reynolds number. Int. J. Multiphase Flow 1992, 18, 145–147. (43) Lun, C. K. K.; Liu, H. S. Numerical simulation of dilute turbulent gas
solid flows in horizontal channels. Int. J. Multiphase Flow 1997, 23 (3), 575–605. (44) Lun, C. K. K. Numerical simulation of dilute turbulent gas
solid flows. Int. J. Multiphase Flow 2000, 26 (10), 1707–1736. (45) Choi, M.; Meisen, A. Hydrodynamics of shallow, conical spouted beds. Can. J. Chem. Eng. 1992, 70 (5), 916–924. (46) Miller, T. F.; Schmidt, F. W. Use of a pressure-weighted interpolation method for the solution of the incompressible Navier
Stokes equations on a nonstaggered grid system. Numer. Heat Transfer, Part A 1988, 14 (2), 213–233. (47) Amsden, A.; Harlow, F. A simplified MAC technique for incompressible fluid flow calculations. J. Comput. Phys. 1970, 6 (2), 322–325. (48) Mathur, K. B.; Epstein, N. Spouted Beds; Academic Press: New York, 1974. (49) Van Velzen, D.; Flamm, H. J.; Langenkamp, H.; Casile, A. Motion of solids in spouted beds. Can. J. Chem. Eng. 1974, 52 (2), 156–161. (50) Olazar, M.; San Jose, M. J.; Alvarez, S.; Morales, A.; Bilbao, J. Measurement of particle velocities in conical spouted beds using an optical fiber probe. Ind. Eng. Chem. Res. 1998, 37 (11), 4520–4527. (51) Olazar, M.; San, J. M. J.; Izquierdo, M. A.; de Salazar, A. O.; Bilbao, J. Effect of operating conditions on solid velocity in the spout, annulus and fountain of spouted beds. Chem. Eng. Sci. 2001, 56 (11), 3585–3594. 4105

dx.doi.org/10.1021/ef200808v |Energy Fuels 2011, 25, 4095–4105