Bubbling Fluidized Bed - American Chemical Society

Apr 3, 2014 - ABSTRACT: Particle-scale study of the gas−solid hydrodynamics and the tube erosion in a 3-D bubbling fluidized bed with a tube bundle ...
2 downloads 0 Views 5MB Size
Article pubs.acs.org/IECR

Particle-Scale Investigation of the Hydrodynamics and Tube Erosion Property in a Three-Dimensional (3-D) Bubbling Fluidized Bed with Immersed Tubes Shiliang Yang, Kun Luo, Jianren Fan,* and Kefa Cen State Key Laboratory of Clean Energy Utilization, Zhejiang University, Hangzhou 310027, People’s Republic of China ABSTRACT: Particle-scale study of the gas−solid hydrodynamics and the tube erosion in a 3-D bubbling fluidized bed with a tube bundle is conducted in the framework of computational fluid dynamics coupled with the discrete element method. General gas−solid flow behaviors are investigated and the time-averaged properties are obtained to explore the influence of tube configuration on the bed hydrodynamics. Moreover, the erosion patterns of immersed tubes locating in the different regions of system are analyzed. The effect of tube configuration on the erosion property is considered. The results show that the presence of the tube bundle limits the bubble diameter and results in homogeneous distribution of the bubble phase in the tube region. Moreover, inserting the tube bundle into the bed reduces the vertical velocity of the solid phase but results in a uniformly horizontal distribution of its concentration. The vertical intensity of solid dispersion is several times that of the lateral one, and the insertion of tube bundle enhances the solid mixing procedure and enlarges the solid dispersion intensities in both the vertical and lateral directions. Different erosion patterns appear on the immersed tubes locating in different regions of the bed. The intrinsic mechanism behind the erosion distribution can be identified from the solid flux distribution near the tubes. using the approach of computational fluid dynamics coupled with the discrete element method (CFD-DEM).14−16 Detailed information on emulsion phase can be captured, giving rise to the ability of deeply exploring the mechanisms governing the solid motion and erosion property in the fluidizing apparatus at the particle-scale level. In this coupling framework, Rong et al.17 conducted a numerical study to analyze the effects of superficial velocity and tube arrangement on the bed hydrodynamics in the fluidized bed with immersed tubes. They found that the impacts between particles and tubes mainly concentrate in the initial period of the bubble wake attack. In the situations with elevated pressure and temperature, Rong and Horio18 numerically investigated the behaviors of particles and bubbles around the immersed tubes in a 2-D fluidized bed using the CFD-DEM coupling approach. The results demonstrated that the bubble size is reduced while the bubble number increases particularly in the tube bank region at the elevated pressure. Gui et al.19 presented a numerical simulation based on the Eulerian− Lagrangian method to investigate the gas−solid hydrodynamics and the particle-tube impact frequency in the bubbling fluidized bed with two immersed tubes. However, only information on the collision frequency on the tube surface is not sufficient to explore the erosion distribution property around the circumferential position of immersed tube. As has been pointed out by Finnie,20 the erosion quantity removed by a colliding particle depends on the particle velocity and the interaction angle. Thus, in order to understand the intrinsic mechanism behind the erosion behavior, there is a need to quantitatively describe the erosion taking into account of the impact angle between the

1. INTRODUCTION Fluidized beds have been widely used in many industrial processes, such as oil refinery, coating, coal gasification/ combustion, and other applications. In order to obtain the desired system performance, a tube bundle is usually introduced into the apparatus.1,2 The presence of a tube bundle alters the gas−solid motion in the system. Simultaneously, the tube surface sustains severe erosion because of intensive collision by the moving particles, which shortens the lifetime of the tube bundle and affects the heat-transfer coefficient of the system. Thus, there is a need of deeply exploring the influence of tube bundle on the bed hydrodynamics and its erosion property for the purpose of achieving an optimal performance. In the past years, many experimental efforts have been taken to investigate the several important aspects of fluidized bed with immersed tubes, such as the heat transfer property of tube surface,3 the bubble hydrodynamics,4 tube erosion,5−7 and the bed expansion behavior.8 In the last two decades, because of the rapid growth in the computational capacity, numerically modeling the gas−solid motion has become a powerful tool to study the bed hydrodynamics in the bubbling fluidized bed.9−11 In general, two methods are widely adopted to track the dense two-phase flow, namely, the Eulerian−Eulerian approach and the Eulerian−Lagrangian approach. In the Eulerian−Eulerian method, both the bubble and emulsion phases are considered as the continuous phases, for which the time averaged Navier− Stokes equations are used to describe their motions. By means of this approach, numerical investigations on the gas−solid hydrodynamics and the tube erosion in the bubbling fluidized bed have been reported.1,2,12,13 Different from the former approach, in the Eulerian− Lagrangian approach, the emulsion phase is tracked at the particle-scale level, and the gas−solid motions are resolved © 2014 American Chemical Society

Received: Revised: Accepted: Published: 6896

September 14, 2013 March 12, 2014 April 3, 2014 April 3, 2014 dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

particle-scale level has been rarely reported in the literature. Besides, to the authors’ knowledge, there has been no any experimental or numerical report on the influence of tube bundle on the solid dispersion behavior in the bubbling fluidized bed. Thus, the current work investigates the hydrodynamics and the erosion distribution around the tube surface in a 3-D bubbling fluidized bed with dense tube bundle by means of the CFD-DEM coupling approach. The gas motion is solved at the computational grid level while the solid motion is tracked individually using the discrete element method. Furthermore, the Finnie erosion model20 is adopted to predict the erosion quantity around the tube surface. General flow behaviors of gas−solid phases are initially explored, followed by the identification of preferential bubble path. Then, quantitative descriptions of the gas−solid hydrodynamics are presented. Finally, the erosion distributions around each tube in the systems with staggered and inline tube configurations are discussed.

Table 1. Equations of Calculating the Stiffness and Damping Coefficients in the Soft Sphere Contact Model Normal stiffness coefficient (kn,ij) and damping coefficient (γn,ij)

kn , ij =

4 Y * R *δn , ij , 3

γn , j = 2

5 β Sn , ijm* 6

Tangential stiffness coefficient (kt,ij) and damping coefficient (γt,ij)

kt , ij = 8G* R *δn , ij ,

γt , ij = 2

5 β St , ijm* 6

in which

Sn , ij = 2Y * R *δn , ij ,

St , ij = 8G* R *δn , ij ,

(1 − νj 2) (1 − νi ) 1 = + , Y* Yi Yj 2

1 1 1 = + , R* Ri Rj

1 1 1 = + m* mi mj β=

ln(e) 2

ln (e) + π

2

,

2(2 + νj)(1 − νj) 2(2 + νi)(1 − νi) 1 = − Yi Yj G*

2. GOVERNING EQUATIONS 2.1. Governing Equations for Fluid Motion. The fluid motion is resolved at the computational grid level and governed by the Navier−Stokes equations. In the framework of large eddy simulation, the mass and momentum equations taking the

colliding particle and the tube surface. Furthermore, exploring the influence of tube configuration on the tube erosion at the

Figure 1. Geometrical description and grid representation of the calculation domain (values given in panels a and b are expressed in units of mm): (a) sketch map of the bed with an inline tube bundle; (b) sketch map of the bed with a staggered tube bundle; (c) grid representation around an inline tube bundle; and (d) grid representation around a staggered tube bundle. 6897

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

Figure 2. Comparison between the calculated results and the experimental data in the literature:31,32 (a) the pressure drop profile with increasing or decreasing the superficial velocity; (b) profile comparison of the time-averaged voidage at height of 0.0072 m, Uf = 0.6 m/s; (c) profile comparison of the timeaveraged voidage at height of 0.0164 m, Uf =0.6 m/s; and (d) profile comparison of the time-averaged solid axial velocity at height of 0.02 m, Uf = 0.6 m/s.

In the momentum equation (eq 2), τ̃f,ij is the filtered stress tensor, which is formulated as

voidage and the interphase momentum exchange into account are formulated as19,21 ∂(ρf εf ) ∂t

+

∂(ρf εf uf̃ , i) ∂t

+

∂(ρf εf uf̃ , i) ∂xi

=0

⎛ ∂uf̃ , i ∂uf̃ , j ⎞ ∂u ̃ ⎟⎟ − 2 μf f , k δij + τf̃ , ij = μf ⎜⎜ 3 ∂xk ∂xi ⎠ ⎝ ∂xj

(1)

∂(ρf εf uf̃ , iuf̃ , j)

= ρf εf g − εf

where μf is the fluid viscosity coefficient and δij is the Kronecker function.

∂xj ∂pf̃ ∂xi

+

∂(εf (τf̃ , ij + τf̃ , ij SGS)) ∂xj

n



∑ i=1

Fd , i ΔV

͠

τ fSGS ̃ , ij (=ρf (uf̃ , iuf̃ , j − uf , iuf , j)) is the subgrid stress tensor, which is modeled based on the eddy viscosity concept. For its calculation, the famous Smagorinsky model22 is adopted in the current work. It can be expressed as

(2)

where ρf, ũf,i, and p̃f are, respectively, the density, the velocity component, and the pressure of fluid phase; n is the total number of particles locating in the current cell; t is the time instant; g is the gravitational acceleration; Fd,i is the drag force exerted on particle i and will be discussed later; and ΔV is the volume of current computational cell. εf is the voidage of fluid phase in the current cell and can be calculated as

εf = 1 −

(4)

⎛ ⎞ 1 τ fSGS ̃ , ij = ρf ⎜2νf , t Sf̃ , ij + τf̃ , kkδij⎟ ⎝ ⎠ 3

̃ is the deformation tensor of the filtered field (S̃f,ij = where Sf,ij (1/2)(∂ũf,i/∂xj + ∂ũf,j/∂xi)). vf,t is the eddy viscosity coefficient, which is calculated as

n ∑i = 1 Vpi

ΔV

(5)

(3)

νf , t = (CsΔ)2 (2Sf̃ , ijSf̃ , ij)1/2

where Vpi is the volume of particle i in the current cell. 6898

(6)

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

Figure 3. Snapshots of bubble motion in the time range of t = 6.0−6.06 s (the bubble boundary is identified with a threshold value of voidage equal to 0.75), Uf = 0.15 m/s: (a−d) system without a tube bundle; (e−h) system with an inline tube bundle.

where Δ is the characteristic length scale used to filter variables (Δ = (Δx × Δy × Δz)1/3; here, Δx, Δy and Δz are the mesh sizes of the computational cell in the X-, Y-, and Z-directions, respectively. Cs is a constant in the Smagorinsky model. It is evaluated using the correlation proposed by Lily,23 as

Cs =

1⎛ 2 ⎞ ⎜ ⎟ π ⎝ 3Ck ⎠

tracked at the particle-scale level. The motion of each particle is governed by Newton’s second law. The equations describing the translational and rotational movements of a particle are formulated as25 mi

3/4

dvi = Fd , i + Fp , i − dt

k

∑ (Fcn, ij + Fdn, ij + Fct , ij + Fdt ,ij) + mi g j=1

(8)

(7)

where Ck is the Kolmogorov constant (Ck = 1.4). In order to capture the eddy in the near-wall region, a wall function of three-layer log-law24 is adopted in the current research. 2.2. Governing Equations for Solid Motion. In the framework of the discrete element method, solid motion is

Ii

dωi = dt

k

∑ (R i × (Fct ,ij + Fdt ,ij) − μr ,ij j=1

× |Fcn , ij|ω̂i

)

(9)

where mi is the mass of particle i, vi the translational velocity of particle i, ωi the rotational velocity of particle i, and Ii the moment of inertia of particle i. The parameter k is the total 6899

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

Figure 4. Snapshots of solid motion in the time range of t = 6.0−6.06 s, Uf = 0.15 m/s: (a−d) system without a tube bundle; (e−h) system with an inline tube bundle.

and wall can be treated in the same manner as dealing with particle−particle collisions. 2.3. Submodels for the CFD-DEM Coupling. 2.3.1. Collision Dynamics Model. In the dense gas−solid motion, a particle frequently collides with other particles or walls. In order to correctly mimic the colliding procedure, the soft sphere contact model27 is chosen in the current work due to its ability to deal with multiple collisions at the same time. The calculation of interaction force can be formulated as

number of particles and walls colliding with the current one. Fcn,ij, Fdn,ij, Fct,ij, and Fdn,ij are, respectively, the normal contact force, the normal damping force, the tangential contact force and the tangential damping force, respectively. Their calculation details will be described in the next section, Submodels for the CFD-DEM Coupling. The parameter μr,ij represents the rolling friction coefficient and ω̂ i stands for the unit angular velocity (ω̂ i = ωi/|ωi|). Ri is the vector from the center of mass of particle i to the contact point of particle i and particle j. Fp,i is the pressure gradient force (Fp,i = −Vp,i∇p̃f),26 since the pressure difference across a surface implies a force difference in the direction of pressure gradient. When a particle collides with the wall, the wall is treated as a stagnant particle of infinite radius and mass. Thus, the collision between a particle 6900

Fcn , ij = kn , ijδn , ij n

(10)

Fdn , ij = γn , ij vn , ij

(11)

Fct , ij = kt , ijδt , ij t

(12)

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

Table 2. Physical and Numerical Parameters Adopted in the Simulation parameter

value Gas Phase

temperature viscosity superficial gas velocity molecular weight pressure density

298 K 1.8 × 10−5 Pa s 0.15 m/s 28.8 kg/mol 1 atm 1.205 kg/m3 Particles

number of particles bed thickness of 1.38 mm bed thickness of 2.76 mm diameter density Poisson ratio Young’s modulus friction coefficient restitution coefficient initial temperature rolling friction coefficient Geometry width thickness height number of cells in the X-direction number of cells in the Y-direction number of cells in the Z-direction total number of grid cells without a tube bundle staggered tube bundle inline tube bundle tube diameter tube temperature tube Vickers hardness

68 696 137 392 0.23 mm 2700 kg/m3 0.25 5.5 × 108 Pa 0.3 0.97 298 K 0.01 20 mm 1.38 mm 100 mm 40 6 130 31 200 30 432 30 336 2.6 mm 298 K 294 MPa

Figure 5. Time statistical property of the solid phase in the system with an inline tube configuration: (a) time evolutionary profile of the average vertical velocity of all the particles in the system, t = 0−7.4 s; (b) profile comparison of time-averaged concentration of solid phase with different statistical time intervals, Z = 0.007 m.

and tangential directions. n and t represent the unit vectors in the normal and tangential directions between the contacting particles, respectively. vn,ij and vt,ij stand for the normal and tangential components of the relative velocity between the colliding particle pair, respectively. The parameters kn,ij and kt,ij represent the stiffness coefficients in the normal and tangential directions, respectively. γn,ij and γt,ij are the damping coefficients in the normal and tangential directions, respectively. The calculation details of stiffness and damping coefficients are listed in Table 1. 2.3.2. Drag Force Model. The drag force exerted by the fluid phase on a suspending particle should take into account of the presence of neighboring particles. It can be formulated as25

Table 3. Parameters Settings in the Bubbling Fluidized Bed Used for Model Validation parameter Domain Size width−depth−height, m grid (width−depth−height) static bed height, m Solid Properties particle number particle density, kg/m3 particle radius, mm Young’s modulus, MPa Poisson ratio restitution coefficient friction coefficient Gas Properties gas density, kg/m3 gas viscosity, kg/(m s)

Fdt , ij = γt , ij vt , ij

value 0.044−0.01−0.2 15−6−68 0.03 9240 1000 1.2 0.12 0.33 0.97 0.1

Fd =

Vpβgs εp

(u f − vp)

(14)

where εp is the solid concentration in the current cell (εp = 1 − εf). uf is the velocity vector of fluid phase in the current cell. βgs is the fluid−particle interaction coefficient, for which the correlation proposed by Koch and Hill28 is adopted and can be expressed as

1.225 1.8 × 10−5

(13)

where δn,ij and δt,ij stand for the relative displacements in the normal and tangential directions, respectively. The subscripts n and t are, respectively, the variable components in the normal

βgs = 6901

18μf εf 2εp ⎛ dp2





F0(εp) +

⎞ 1 F3(εp)Rep⎟ ⎠ 2

(15)

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

Figure 6. Influence of bed thickness on the time-averaged property of solid phase, Z = 0.007 m: (a) solid vertical velocity and (b) solid concentration.

where dp is the particle diameter. The Reynolds number (Rep) is calculated as Rep =

εf ρf |u f − vp|dp μf

(16)

Figure 7. Contour plots of time-averaged voidage and its standard variance in the central slice of bed, Y = 0.69 mm, Uf = 0.15 m/s: (a) voidage in a system free of tube bundles; (b) voidage in a system with an inline tube configuration; (c) standard variation of voidage in a system free of tube bundles; and (d) standard variation of voidage in a system with an inline tube configuration.

F0 and F3 are the model parameters used in the model, which can be calculated as28 ⎧ εp 135 ⎪ 1 + 3 2 + 64 εp ln(εp) + 16.14εp εp < 0.4 ⎪ ⎪ 1 + 0.681εp − 8.48εp 2 + 8.16εp3 F0(εp) = ⎨ ⎪ 10εp ⎪ εp ≥ 0.4 ⎪ 3 ⎩ εf

F3(εp) = 0.0673 + 0.212εp +

0.0232 εf 5

where PH is the Vickers hardness of the tube surface. f(γ) is a dimensionless function of the impact angle γ, which is formulated as

(17)

2 ⎧ ⎪ sin(2γ ) − 3 sin (γ ) tan(γ ) < 1/3 f (γ ) = ⎨ ⎪ 2 tan(γ ) > 1/3 ⎩ cos (γ )/3

(18)

2.3.3. Erosion Model. When the immersed tube is exposed to the dense two-phase flow, erosion appears at the tube surface, because of the vigorous collisions by the moving particles. The erosive quantity by a moving particle mainly depends on the impact velocity and the impact angle of the collision. In this work, the Finnie erosion model20 is adopted to predict the tube erosion. It is formulated as 1 E= mp vp 2f (γ ) 8PH (19)

(20)

3. COUPLING PROCEDURE AND SIMULATION CONDITIONS 3.1. The Coupling Procedure between CFD and DEM. The CFD-DEM coupling approach has been well-documented in the literature.25 The current work extends the coupling procedure with erosion prediction. The coupling occurs between the CFD solver and DEM solver, which are used to solve the governing equations of fluid motion and solid motion, 6902

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

Figure 8. Profile comparison of the time-averaged solid vertical velocity in the systems with different tube configurations, Y = 0.0069 m: (a) Z = 0.007 m; (b) Z = 0.017 m; (c) Z = 0.022 m; and (d) Z = 0.027 m.

Figure 9. Profile comparison of the time-averaged concentration of solid phase in the systems with different tube configurations, Y = 0.0069 m: (a) Z = 0.007 m, (b) Z = 0.017 m, (c) Z = 0.022 m, and (d) Z = 0.027 m.

momentum exchange combined with the voidage obtained. Then, the governing equations of fluid phase are solved. Since the velocity and pressure couple with each other through the momentum equation of fluid motion, the iteration of governing

respectively. At each time step, position information on each particle is used to detect the computational cell that this particle locates in, and then voidage of each computational cell is obtained. Solid velocity is used to estimate the interphase 6903

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

Figure 10. Qualitative comparison of the mixing process in systems with different tube configurations.

equations is carried out using the Pressure Implicit with Splitting of Operator (PISO) algorithm,29 which consists of one predictor step and several corrected steps. After this, DEM solver estimates the interaction forces between the fluid and solid phases, and then the collision detection is applied to calculate the collision forces in the colliding particle−particle or particle−wall pair. The velocity and position of a particle are calculated by integrating the Newtown’s second law. In the tracking procedure, erosion calculation is triggered if the collision between a particle and tube surface occurs. Then, DEM solver prepares the information on solid phase for next calculation time step. 3.2. Simulation Conditions. The investigated bubbling fluidized bed has a section of 20 mm × 1.38 mm (width × depth) and a height of 100 mm. In order to investigate the influence of the tube bundle on the gas−solid motion, both staggered and inline tube bundles are inserted into the system, respectively. Each tube has a diameter of 2.6 mm in the calculation and is identified with a specific number. The geometrical description of the calculation domain with the tube bundle is schematically illustrated in Figure 1. The system without the tube bundle is divided into 40 × 6 × 130 elements in the X-, Y-, and Z-directions, respectively. An unstructured grid is adopted to capture the tube surface. The total cell numbers for the systems with staggered and inline tube bundles are 30 432 and 30 336, respectively. Furthermore, the circumferential positions of each tube are described in a clockwise

direction with 0° at the top of the tube. Detailed physical and numerical parameters used in the current work are listed in Table 2. In order to solve the governing equations of fluid motion, unsteady-state items in equations are discretized with the Crank− Nicholson scheme, and the diffusion items are discretized with the Gauss linear scheme, while the first-order upwind scheme is used for the convective items. For the DEM solver of the solid phase, the position and velocity of each individual particle are obtained with explicitly integrating the governing equations. 3.3. Boundary and Initial Conditions. Different boundary conditions are used for variables in different regions of the bed. In order to reduce the computational capacity, the front and back of the system are assigned with the periodic boundary conditions for all the variables involved in both the CFD and DEM aspects. For the velocity, the uniform distribution is applied to the inlet of fluidized bed; the Neumann boundary conditions are assigned to the outlet. For the boundary condition of pressure, the environmental atmospheric pressure is adopted for the outlet, while other boundaries are assigned the Neumann boundary conditions. For the initial distribution of particles in the system, a total number of 68 696 particles are scattered randomly in the bed, and then fall under the influence of gravity in the presence of an immersed tube. After the total kinematic energy of the entire system decays to nearly zero, the fluidizing gas is introduced from the bed inlet. In order to make the 6904

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

Figure 11. Quantitative evaluation of mixing intensity in the systems with different tube configurations: (a) illustration of sampling cells adopted to calculate the mixing index; (b) time evolutionary profiles of the mixing index in the systems with different tube configurations.

model capture the important dynamical property of the bubbling fluidized bed. Based on the proposed model, the hydrodynamics and solid dispersion in the bubbling fluidized bed are investigated.33,34 Moreover, it has been adopted to explore the solid mixing and process identification in the internally circulating fluidized bed.11,16,21 4.2. General Flow Behavior of Gas Phase. The influence of the tube bundle on the gas−solid motion is qualitatively obtained by comparing the related aspects in the system free of a tube bundle and the system with an inline tube bundle. Figure 3 shows snapshots of bubble motion in the time range of 6.0−6.06 s. The boundary of the bubble phase is identified with a threshold value of voidage equal to 0.75. In the bed without a tube bundle, small bubbles are generated randomly at the bed bottom and then rise upward. In the rising processes, interactions between neighboring bubbles appear. Bubble coalescence in both the vertical and lateral directions frequently appears due to the instability of bubble boundary, giving rise to the bubbles with large diameter in the upper region of system. Comparatively, the tube bundle explicitly affects the general flow behavior of gas phase. In the system with inline tube configuration, the bubble grows until the first row of tube bundle is reached. Then, some bubbles are split into two or even more smaller bubbles. For the large bubble, the breakup will not occur and the tube is wholly surrounded by the bubbles, followed by the bubble being elongated. After passing the lowest row of the tube bundle, bubble coalescence appears again before the next row is encountered. This procedure repeats until the region free of tubes is reached. Thus, the presence of tube bundle suppresses the growth of the bubble, leading to a limited bubble size in the whole bed. Small bubbles exist in the system at the same time with a relatively uniform

calculation stable, the time step of calculation should be small enough. In order to accurately mimic the collision procedure, the time step for the solid motion should be smaller than the critical value.30 In the current work, the time step for fluid motion is 1 × 10−5 s, while it is 1 × 10−7 s for solid motion. The coupling between CFD and DEM occurs at every 100 time steps of solid motion. In order to construct the collision list for a particle, the detection distance is chosen to be twice the particle diameter in the calculation. Each simulation case is conducted for a real time of 15 s.

4. RESULTS AND DISCUSSION 4.1. Model Validation. Validation of the proposed model is carried out based on the experiment conducted by Müller et al.31,32 The investigated bubbling fluidized bed has dimensions of 44 mm × 10 mm × 200 mm along the width, depth, and height, respectively. The formed packed bed with a height of 0.03 m consists of 9240 particles. The validation is mainly achieved by comparing the minimum fluidization velocity (Umf) and the time-averaged bed hydrodynamics. Details of the parameters adopted in the simulation are listed in Table 3. In order to test the minimum fluidization velocity, the pressure drop of the system with gradually increasing or decreasing the superficial velocity is recorded. The minimum fluidizing velocity can be identified with the critical value corresponding to the sudden change of pressure drop. As illustrated in Figure 2a, the predicted value of this value is 0.3 m/s, which is the same as that obtained in the experiment.31 Furthermore, the timeaveraged value of bed hydrodynamics are obtained in the system fluidized at the superficial velocity of 0.6 m/s (Uf = 2Umf). As illustrated in Figures 2b−d, the predicted profiles of timeaveraged axial velocity of solid phase and voidage of gas phase agree well with the experimental data, indicating the proposed 6905

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

distribution across the bed. Furthermore, it can be obtained that there is a bubble jet existing in the space between the tubes, because of the area reduction of the rising tunnel, which has been also reported in the experimental work of Wiman et al.,6 as well as in the numerical work of Costa et al.35 4.3. General Flow Behavior of Solid Phase. Particles in the fluidized bed undergo vigorous motion with rising bubbles. Figure 4 shows the snapshots of solid motion in the time range of 6.0−6.06 s. Some particles are entrained into the wake of the formed bubble in the bed bottom and then rise upward. Vigorous upward motion of solid phase mainly appears in the central region of the bed especially in the dome and wake of bubble phase. After the eruption of bubble phase at the top region, particles are scattered into the free domain of the bed, giving rise to the intensive downward flow near the bed wall. In the system with the tube bundle, when the tube is encountered, particles in the bubble nose move in an opposite direction after colliding the tube surface and then rain down through the bubble, giving rise to intensive collision between the falling particles and the upward particles in the bubble wake. The others mainly concentrate in the rising tunnel formed in the gaps between the inline tube columns, resulting from small flow resistance in the vertical direction. It is clearly observed that the immersed tubes are always surrounded by particles. The occasional formation of the stagnant cap in the top region of each tube is frequently swept by the rising bubbles. This trend agrees well with the previous experimental investigations36 and numerical results.13,17,37 4.4. Time Statistical Properties. 4.4.1. Influence of Statistical Time and Bed Thickness. In order to quantitatively investigate the bed hydrodynamics, the time-averaged procedure should be carried out with enough time. Besides, since the periodic boundary condition is assigned for the front and back of the investigated geometry, there is a need to explore the sufficiency of bed thickness adopted to apply this boundary condition. Based on the system with inline tube configuration, these two aspects are investigated. After introducing the gas phase, particles in the system are fluidized. A startup procedure appears until the pseudo-steady-state operation is reached. The time evolutionary profile of the averaged vertical velocity of all the particles in the system is obtained. As illustrated in Figure 5a, the profile fluctuates due to the turbulent and chaotic motion of solid phase in the fluidizing procedure. At the initial stage of calculation, large fluctuation of solid vertical velocity can be obtained, which corresponds to the startup procedure. After a small period of adjustment, the investigated variable seems to reach the pseudo-steady-state status. Conservatively, the time-averaged procedure is carried out from the time instant of t = 5 s. Furthermore, it is carried out with three different time intervals. As shown in Figure 5b, the profile becomes smooth with a total statistical time of 10 s. Thus, to investigate the time-averaged property in the current work, the data of initial 5 s are discarded and time-averaged analysis is conducted over the data of last 10 s in all the simulated cases. In order to evaluate the sufficiency of bed thickness for utilizing the periodic boundary condition, the time-averaged concentration and vertical velocity of solid phase in the system with two different bed thicknesses are compared and shown in Figure 6. General tendency of these profiles behaves in a similar pattern and minor difference appears between them, indicating that the adopted thickness is proper.

Figure 12. Time evolutionary profiles of (a) vertical and (b) lateral dispersion coefficients of the solid phase in the bed free of tube bundles, t = 5−15 s.

4.4.2. Preferential Bubble Path. Figure 7 illustrates the contour plots of the average and standard variance of voidage in the central slice (Y = 0.69 mm) of the systems with or without immersed tubes. It can be observed that large voidage mainly concentrates in the central region of the bed, because of the restriction of the wall effect on gas motion. A relatively small amount of voidage exists in the near-wall region, because of the return flow generated by the falling particles after the bubble eruption. Furthermore, the insertion of tube bundle results in the expansion of bed height. It has been mentioned in the literature5,6,38 that there is a preferential bubble path that is mainly focusing on the region with large values of both the average and standard variance of voidage. In the system without a tube bundle, the preferential path of bubble motion mainly lies in the central region of the bed. Moreover, the flow area of bubble phase decreases along the bed height. The insertion of a tube bundle results in a more homogeneous distribution of the bubble phase in the region with the tube bundle. It should be noted that the large variance of voidage in the regions below tubes is mainly due to the backward movement of rising particles after colliding with the internal surface, indicating that particles frequently migrate near the bottom region of the tube. 4.4.3. Statistic Properties of Solid Phase. In order to study the influence of tube bundle on the solid motion in the bed, the transversal distributions of solid vertical velocity and concentration in different regions of the bed are compared. The data 6906

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

located in the central slice (Y = 0.69 mm) of the bed. Figure 8 shows the time-averaged profiles of the solid vertical velocity at four investigated heights. In general, the vertical motion of the solid phase mainly concentrates in the central region of the bed, because of the vigorous motion of bubble phase, while strong back-mixing of solid motion appears in the near wall region. At a height of Z = 0.007 m (below the tube bundle), only a small difference between the systems with different tube configurations can be obtained. Compared with the profiles in the bed without a tube, the insertion of the tube bundle limits the vertical motion of the solid phase in the tube region, giving rise to the appearance of several peaks and valleys in the distribution profiles. Compared to the system with a staggered tube configuration, the solid velocity under each individual tube of the system with inline configuration is smaller

Table 4. Time-Averaged Value and Standard Variance of Solid Dispersion Coefficient in the Lateral or Vertical Direction Lateral Dispersion Vertical Dispersion Coefficient (× 10−5 m2/s) Coefficient (× 10−5 m2/s)

free of tube inline tube bundle staggered tube bundle

mean

standard variation

mean

standard variation

1.85 3.50 3.40

0.64 0.89 0.90

5.71 10.44 10.12

1.69 2.73 3.05

are obtained at the heights of Z = 0.007 m (below the tube bundle), Z = 0.017 m (between the first and second row of tubes from bottom to top), Z = 0.022 m (in the middle place of second and third row of tubes), Z = 0.027 m (above the tube bundle),

Figure 13. Erosion distribution around the circumferential positions of the immersed tubes in the inline system (the value stands for the accumulated erosion quantity in the time range of t = 0−15 s), Uf = 0.15 m/s: (a) tube 7; (b) tube 8; (c) tube 4; (d) tube 5; (e) tube 1; and (f) tube 2. 6907

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

initial mixing process, the slope of mixing profile is large, indicating that the mixing intensity is the strongest. After the dynamic mixing equilibrium is established, the mixing index fluctuates around a constant value, which corresponds to the micromixing process. The time required for the bed with tube bundle to reach the equilibrium is 1.25 s, while it is nearly 3.0 s for the bed without a tube bundle. However, the effect of tube configuration on the mixing intensity is not obvious in the current work, mainly because of the small domain that is investigated. 4.6. Effect of Tube Bundle on Solid Dispersion. Solid dispersion in the bubbling fluidized bed is crucial for gas−solid contacting efficiency, especially in the combustion procedure.40 There has been a continuous interest of evaluating the solid dispersion behavior in the fluidized bed using both experimental and numerical approaches.14,34,41,42 However, the influence of tube bundle on the solid dispersion has not been reported. The solid dispersion intensity can be described using the solid dispersion coefficient. Two approaches can be used to evaluate the solid dispersion coefficient, namely the macro method and the micro method.41 The former is based on fitting the particle concentration by using a Fickian-type diffusion equation, while the latter evaluate the dispersion intensity from the particle displacement. By means of a CFD-DEM coupling approach, the latter approach is adopted in the current work. For a specific particle i, the dispersion coefficient in a specific direction can be evaluated, using Einstein’s equations,41 as

while a rising tunnel with a relatively large value of solid vertical velocity is formed in the gap between the tubes. The value difference between the peaks and valleys of the solid vertical velocity is larger for the inline tube arrangement, compared to the staggered one. Figure 9 compares the transversal distributions of solid concentration in the three investigated systems. In general, the solid concentration is small in the central region, because of the large gap, while it is small near the wall, because of the downward motion of the solid phase. At all of these investigated heights, the insertion of a tube bundle results in the uniform distribution of the solid phase across the bed width. Furthermore, the peak of solid concentration appears below the tube bundle, because of the intensive collision of rising particles on the tube surface. When operating under the same conditions, a more uniform distribution of solid concentration can be obtained in the system with staggered tube configuration, compared with the inline one. The differences in the flow hydrodynamics between the systems with different tube configurations are consistent with the experimental report of Wiman et al.,6 as well as the numerical work of Rong et al.,17 Schreiber et al.,2 and Li et al.13 4.5. Effect of Tube Bundle on Solid Mixing. Solid mixing is commonly encountered in the natural and industrial processes. Quantitatively evaluating the tube bundle on solid mixing is extremely useful for the design and operation of the apparatus. Figure 10 illustrates the mixing process of the solid phase in the systems with different tube configurations. In order to intuitively describe the mixing process of the solid phase, the particles in the system are vertically divided into two groups, assigned with different labels before the calculation. Influenced by the vigorous motion of bubble phase in the operating procedure, chaotic motion of the solid phase can be observed. The particles in the wake of bubble phase rise vertically along the bed height while they fall along the wall, leading to the vertical mixing of particles. Furthermore, a part of the rising particles are drifted out from the bubble wake in the rising process, while others are entrained into it, leading to the local mixing of the solid phase. As time passes, the boundary between these two groups diminishes with particles dispersing into each other. Finally, a dynamical mixing procedure is established after sufficient time. As can be vividly obtained by comparing the mixing process in the system with tubes, the insertion of tube bundle into the bed enhances the mixing process, mainly resulting from its ability to break the large cluster consisting of particles from the same group. Quantitative estimation of the mixing intensity in the operation process is achieved by adopting the Lacey mixing index.39 According to the Lacey theory,39 the probability of finding one component of the mixture at any position should be fixed in a completely random mixture, and it should be the same as the overall portion of this component in the mixture. In order to calculate the mixing index, the investigated domain should be divided into plenty of sampling cells. Subsequently, the concentration of a specific particle group in each cell is evaluated, followed by calculating their variance. Details of equations used to calculate the mixing index can be found in our previous work.16 In the current work, as schematically illustrated in Figure 11a, the calculation domain is divided into 10 × 3 × 20 cells in the X-, Y-, and Z-directions, respectively. These sampling cells are orderly labeled but only the cells containing particles are adopted when calculating the mixing index. The time evolutionary profiles of solid mixing in the three investigated systems are shown in Figure 11b. In the

Di =

(Δri)2 2Δt

(21)

where Di is the local dispersion coefficient of particle i along a specific direction. Δri is the local displacement of current particle along the specific direction. Δt is the time interval. The dispersion intensity of the whole system is averaged over all the dispersion coefficients of particles in the system. It is expressed as D=

1 NP

NP

∑ i=1

(Δri)2 2Δt

(i = 1, 2 , ..., NP)

(22)

where NP is the total number of particles in the system. Figure 12 shows the time evolutionary profiles of the lateral and vertical dispersion coefficients of the solid phase in the system without a tube bundle. In general, the dispersion coefficient fluctuates around a constant value due to the chaotic

Figure 14. Vector plot of the time-averaged solid flux in the region with an inline tube bundle (SF stands for the solid flux, measured in units of kg/(m2s)). 6908

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

Figure 15. Erosion distribution around the circumferential positions of the immersed tubes in the staggered system (the value stands for the accumulated erosion quantity in the time range of t = 0−15 s), Uf = 0.15 m/s: (a) tube 6, (b) tube 7, (c) tube 4, (d) tube 1, and (e) tube 2.

and turbulent motion of the solid phase in the bed. Moreover, it can be observed that the vertical dispersion intensity of solid phase is several times greater than the lateral one, mainly due to the vertical introduction of the gas phase into the system. The mean value of the solid dispersion coefficient and its standard variance in the systems with different configurations are presented in Table 4. As can be observed, the dispersion intensity in the system with a tube bundle is nearly twice of that in the system that is free of the tube bundle. Thus, the insertion of the tube bundle enhances the dispersion intensity of solid motion in both the lateral and vertical directions, mainly due to the enhanced chaotic motion of the solid phase in the bed. 4.7. Erosion Distribution around Tubes. The tube bundle is exposed to the dense gas−solid flow, resulting in

erosion at the tube surface. Erosion distribution in the tube bundle is evaluated only for the tubes locating in the left half part of the bed. Figure 13 shows the accumulated erosion quantity around the tube surface in the system with an inline tube configuration. Because of the short simulation time, the fluctuations along the circumference are rather large. Thus, only qualitative interpretation of erosion distribution is carried out. Different erosion patterns appear for the tubes locating in different regions of the system. In general, the erosion mainly takes place at the bottom of each tube, while the top part of the tube suffers little erosion. This is in accordance with the experimental work conducted by Zhu,43 Wiman et al.,5 and Johansson et al.6,7 Besides, the maximum erosion quantity appears at different positions of the tubes. Since the erosion 6909

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

Besides, the insertion of tube bundle enlarges the bed height and results in a homogeneous distribution of bubble phase in the region with the tube bundle. (2) The insertion of tube bundle limits the vertical motion of solid phase in the tube region. As compared with the system of inline tube configuration, the value difference between the peaks and valleys of solid vertical velocity is larger. The presence of tube bundle results in a more uniform distribution of solid concentration across the bed width. (3) The time required to reach the dynamical equilibrium of solid mixing is reduced by mounting the tube bundle in the bed. In the bubbling fluidized bed, the solid dispersion coefficient fluctuates around a constant value. The vertical dispersion intensity of the solid phase is several times greater than the lateral one. Inserting the tube bundle in the system enhances the dispersion intensities of the solid phase in both the vertical and lateral directions. (4) Different erosion patterns appear for the immersed tubes locating in different locations of the system. The key to understand the erosion distribution is the detailed knowledge of the time-averaged solid flux distributed near the tube. In the inline tube configuration, severe erosion appears in the lower region of the lowest tube. For the tubes locating near the wall, large erosion can be observed in the lower part of the tube facing the center of the bed. Except for the tube located in the lowest row of the tube bundle, tubes in the central region sustains severe erosion at the two sides of the lower part. However, the severe erosion appears in the bottom region of the tubes in the system with a staggered tube bundle.

quantity strongly depends on the particle velocity, the impact angle, and solid concentration, these effects can be quantitatively described by the solid flux. The instantaneous solid flux (φ) is obtained by the multiplication of the density ρp, the concentration εp(x, t), and the velocity of solid phase. The time-averaged solid flux is calculated as ⟨φp(x)⟩ =

1 Nt

Nt

∑ (εp(x, t )ρp vp(x, t )) 0

(23)

where x, t, and Nt are the cell number, the time instant, and the total number of time intervals, respectively. Figure 14 presents the vector plot of time-averaged solid flux around the tube surface. At first, the investigation of erosion pattern is mainly focused on the tubes in the near-wall region. Asymmetric erosion pattern appears for these three tubes. The maximum erosion quantity exists on the right side of the lower part of these tubes (at an angle of 120°), mainly due to direct collision by the vigorous moving particles along the rising tunnel.5 However, small erosion exists at the 60° and 240° positions, where large solid flux deviates from the tube surface. The difference appears for distribution of tube 7 locating in the left part of the highest row, as shown in Figure 13a, and the maximum value of erosion quantity appears at the 125° position and while the minimum exists at the 60° position. Significantly large erosion quantity appears in the left upper of this tube, compared with tubes 1 and 4, resulting from the intensive collision by the return flow generated after the bubble eruption at the bed surface. This is in accordance with the experimental report of Wiman et al.5 Directly influenced by the upward motion of solid phase in the central region of the bed, severe erosion appears in the lower part of tube 2. For tubes 5 and 8, located in the central line of the system, the erosion pattern shows a relatively symmetrical distribution with the maximum appears at the 125° and 235° positions of each tube. Different arrangements of the tube bundle alter the gas−solid hydrodynamics in the bed, which can significantly change the erosion patterns of immersed tubes in the system. Figure 15 shows the accumulated erosion quantities for each individual tube in the system with a staggered tube configuration. Except for the tubes in the lowest row, different erosion patterns can be obtained, in comparison with those in the system with an inline configuration. For tubes located in the central line of the bed, the erosion in the lower part of tube 7 shows a more uniform pattern. For tubes located in the second row, the area with large erosion quantity can be observed in the right lower part of the tube (tube 4 in Figure 15). For the tubes located in the left upper region of the bundle (tube 6 in Figure 15), the erosion quantity in the left upper region shows a scale similar to that shown in the right lower part of this tube.



AUTHOR INFORMATION

Corresponding Author

*Fax: +86-0571-87991863. E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS Financial support from the Major Program of the National Natural Science Foundations of China (Grant Nos. 51390491, 51390493), and the National Program for Special Support of Top-Notch Young Professionals are sincerely acknowledged.



5. CONCLUSION In the CFD-DEM coupling framework, numerical simulation of the gas−solid motions in a 3-D bubbling fluidized bed with or without a tube bundle is carried out to analyze the gas−solid hydrodynamics and the tube erosion property. Moreover, the influences of tube configuration on the gas−solid hydrodynamics are discussed. Based on the calculated results, several conclusions can be drawn: (1) The presence of immersed tubes suppresses the bubble growth, leading to a small bubble size in the region of tube bundle. Bubble jets appear in the space between tubes and a stagnant cap appears occasionally at the top of the tube. 6910

NOMENCLATURE Ck = Kolmogorov constant Cs = Smagorinsky constant coefficient D = dispersion coefficient of solid phase, m2/s dp = particle diameter, m e = restitution coefficient E = erosion quantity, m3 F0, F3 = model parameters used in the drag force model Fcn,ij = normal contact force between particle i and particle j, N Fct,ij = tangential contact force between particle i and particle j, N Fd,i = drag force exerted on particle i, N Fdn,ij = normal damping force between particle i and particle j, N Fdt,ij = tangential damping force between particle i and particle j, N Fp,i = far-field pressure force exerted on particle i, N g = gravitational acceleration, m/s2 G* = tangential stiffness coefficient of solid phase dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

γn = damping coefficient in normal direction, kg/s γt = damping coefficient in tangential direction, kg/s τf,ij = viscous stress tensor, Pa τSGS f,ij = subgrid stress tensor, Pa v = Poisson ratio of solid phase, dimensionless

I = particle moment of inertia, kg m2 k = total number of particles in contact with the current one kn,ij = normal stiffness of solid phase, N/m kt,ij = tangential stiffness of solid phase, N/m m = particle mass, kg m* = effective mass of a particle, kg Np = total number of particles in the system Nt = number of time intervals for statistics n = normal unit vector between colliding particles n = number of particle locating in the current cell pf = pressure, Pa PH = Vickers hardness of tube surface, Pa R* = effective radius of a particle, m Ri = vector from the mass center of particle i to the contact point of particle i and particle j, m Rep = particle Reynolds number Δri = particle displacement, m Sij = deformation tensor of the filtered field or resolved strain rate, 1/s Sn,ij = normal damping coefficient St,ij = tangential damping coefficient t = time instant, s t = tangential vector between two colliding particles Δt = time interval, s t = tangential unit vector between colliding particles Uf = fluid velocity, m/s uf = velocity vector of the fluid phase in the computational cell, m/s Umf = minimum fluidization velocity, m/s uf,i, uf,j = velocity components of the fluid phase, m/s Usz = vertical velocity of the solid phase, m/s vi = particle velocity, m/s vn,ij = normal component of the relative velocity between colliding pairs, m/s vt,ij = tangential component of the relative velocity between colliding pairs, m/s Vp = particle volume, m3 ΔV = volume of the current cell, m3 x = cell number X, Y, Z = coordinate axes, dimensionless Yp = actual Young’s modulus Y* = effective Young’s modulus

Subscripts

c = contact force d = drag force f = fluid phase i = particle i j = particle j n = normal component of variable p = particle phase t = tangential component of variable Abbreviations



CFD = computational fluid dynamics DEM = discrete element method PISO = Pressure Implicit with Splitting of Operator SF = solid flux 2-D = two-dimensional 3-D = three-dimensional

REFERENCES

(1) Yusuf, R.; Halvorsen, B.; Melaaen, M. C. Eulerian−Eulerian simulation of heat transfer between a gas−solid fluidized bed and an immersed tube-bank with horizontal tubes. Chem. Eng. Sci. 2011, 66 (8), 1550−1564. (2) Schreiber, M.; Asegehegn, T. W.; Krautz, H. J. Numerical and Experimental Investigation of Bubbling Gas−Solid Fluidized Beds with Dense Immersed Tube Bundles. Ind. Eng. Chem. Res. 2011, 50 (12), 7653−7666. (3) Masoumifard, N.; Mostoufi, N.; Hamidi, A.; Sotudeh-Gharebagh, R. Investigation of heat transfer between a horizontal tube and gas− solid fluidized bed. Int. J. Heat Fluid Flow 2008, 29 (5), 1504−1511. (4) Lyczkowski, R. W.; Bouillard, J. X.; Gamwo, I. K.; Torpey, M. R.; Montrone, E. D. Experimental and CFD Analyses of Bubble Parameters in a Variable-Thickness Fluidized Bed. Ind. Eng. Chem. Res. 2009, 49 (11), 5166−5173. (5) Wiman, J.; Mahpour, B.; Almstedt, A. E. Erosion of horizontal tubes in a pressurized fluidized bedInfluence of pressure, fluidization velocity and tube-bank geometry. Chem. Eng. Sci. 1995, 50 (21), 3345−3356. (6) Wiman, J.; Almstedt, A. E. Hydrodynamics, erosion and heat transfer in a pressurized fluidized bed: Influence of pressure, fluidization velocity, particle size and tube bank geometry. Chem. Eng. Sci. 1997, 52 (16), 2677−2695. (7) Johansson, K.; Norling, R.; Hjörnhede, A.; Almstedt, A.; Johnsson, F.; Nylund, A. Hydrodynamics and steel tube wastage in a fluidized bed at elevated temperature. Chem. Eng. Sci. 2004, 59 (1), 31−40. (8) Löfstrand, H.; Almstedt, A. E.; Andersson, S. Dimensionless expansion model for bubbling fluidized beds with and without internal heat exchanger tubes. Chem. Eng. Sci. 1995, 50 (2), 245−253. (9) Hulme, I.; Clavelle, E.; van der Lee, L.; Kantzas, A. CFD Modeling and Validation of Bubble Properties for a Bubbling Fluidized Bed. Ind. Eng. Chem. Res. 2005, 44 (12), 4254−4266. (10) Li, T.; Mahecha-Botero, A.; Grace, J. R. Computational Fluid Dynamic Investigation of Change of Volumetric Flow in Fluidized-Bed Reactors. Ind. Eng. Chem. Res. 2010, 49 (15), 6780−6789. (11) Fang, M.; Luo, K.; Yang, S.; Zhang, K.; Fan, J. LES−DEM investigation of gas−solid flow dynamics in an internally circulating fluidized bed. Chem. Eng. Sci. 2013, 101, 213−227. (12) Gustavsson, M.; Almstedt, A. E. Two-fluid modelling of coolingtube erosion in a fluidized bed. Chem. Eng. Sci. 2000, 55 (4), 867−879.

Greek Symbols

βgs = interphase momentum transfer coefficient, kg/(m3 s) γ = impact angle between the particle and surface, ° δij = Kronecker delta δn,ij = normal displacements between particle i and particle j, m δt,ij = tangential displacements between particle i and particle j, m Δ = subgrid characteristic length scale, m Δx = mesh size in the X direction, m Δy = mesh size in the Y direction, m Δz = mesh size in the Z direction, m εf = voidage εp = solid concentration μf = gas dynamic viscosity, kg/(m s) μr,ij = rolling friction coefficient ρf = gas density, kg/m3 νt = eddy viscosity coefficient, kg/(m s) φ = solid flux, kg/(m2 s) ω = particle angular velocity, 1/s ω̂ i = unit angular velocity 6911

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912

Industrial & Engineering Chemistry Research

Article

Part II: Solid dispersion and circulation properties. Powder Technol. 2014, 256, 395−403. (35) Costa, A. M. S.; Colman, F. C.; Paraiso, P. R.; Jorge, L. M. M. CFD modelling of fluidized bed with immersed tubes. Adv. Chem. Eng. 2012, 357−378. (36) Wong, Y. S.; Seville, J. P. K. Single-particle motion and heat transfer in fluidized beds. AIChE J. 2006, 52 (12), 4099−4109. (37) Schmidt, A.; Renz, U. Numerical prediction of heat transfer between a bubbling fluidized bed and an immersed tube bundle. Heat Mass Transfer 2005, 41 (3), 257. (38) Li, T.; Dietiker, J.; Zhang, Y.; Shahnam, M. Cartesian grid simulations of bubbling fluidized beds with a horizontal tube bundle. Chem. Eng. Sci. 2011, 66 (23), 6220−6231. (39) Lacey, P. M. C. Developments in the theory of particle mixing. J. Appl. Chem. 1954, 4 (5), 257−268. (40) Pallarès, D.; Johnsson, F. A novel technique for particle tracking in cold 2-dimensional fluidized bedsSimulating fuel dispersion. Chem. Eng. Sci. 2006, 61 (8), 2710−2720. (41) Liu, D.; Chen, X. Lateral solids dispersion coefficient in largescale fluidized beds. Combust. Flame 2010, 157 (11), 2116−2124. (42) Olsson, J.; Pallarès, D.; Johnsson, F. Lateral fuel dispersion in a large-scale bubbling fluidized bed. Chem. Eng. Sci. 2012, 74, 148−159. (43) Zhu, J. Tube Erosion in Fluidized Beds; The University of British Columbia: Vancouver, BC, Canada, 1982.

(13) Li, T.; Dietiker, J.; Zhang, Y.; Shahnam, M. Cartesian grid simulations of bubbling fluidized beds with a horizontal tube bundle. Chem. Eng. Sci. 2011, 66 (23), 6220−6231. (14) Luo, K.; Yang, S.; Zhang, K.; Fang, M.; Fan, J. Particle Dispersion and Circulation Patterns in a 3D Spouted Bed with or without Draft Tube. Ind. Eng. Chem. Res. 2013, 52 (28), 9620−9631. (15) Yang, S.; Luo, K.; Fang, M.; Zhang, K.; Fan, J. ThreeDimensional Modeling of Gas−Solid Motion in a Slot-Rectangular Spouted Bed with the Parallel Framework of the Computational Fluid Dynamics−Discrete Element Method Coupling Approach. Ind. Eng. Chem. Res. 2013, 52 (36), 13222−13231. (16) Fang, M.; Luo, K.; Yang, S.; Zhang, K.; Fan, J. Computational Fluid Dynamics−Discrete Element Method Investigation of Solid Mixing Characteristics in an Internally Circulating Fluidized Bed. Ind. Eng. Chem. Res. 2013, 52 (22), 7556−7568. (17) Rong, D.; Mikami, T.; Horio, M. Particle and bubble movements around tubes immersed in fluidized bedsA numerical study. Chem. Eng. Sci. 1999, 54 (23), 5737−5754. (18) Rong, D.; Horio, M. Behavior of particles and bubbles around immersed tubes in a fluidized bed at high temperature and pressure: A DEM simulation. Int. J. Multiphase Flow 2001, 27 (1), 89−105. (19) Gui, N.; Fan, J. R.; Luo, K. DEM−LES study of 3-D bubbling fluidized bed with immersed tubes. Chem. Eng. Sci. 2008, 63 (14), 3654−3663. (20) Finnie, I. Erosion of surfaces by solid particles. Wear 1960, 3 (2), 87−103. (21) Luo, K.; Fang, M.; Yang, S.; Zhang, K.; Fan, J. LES−DEM investigation of an internally circulating fluidized bed: Effects of gas and solid properties. Chem. Eng. J. 2013, 228, 583−595. (22) Smagorinsky, J. General circulation experiments with the primitive equations. I. The basic experiment. Mon. Weather Rev. 1963, 91 (3), 99−164. (23) Lilly, D. K. A proposed modification of the Germano subgridscale closure method. Phys. Fluids 1992, 4, 633−635. (24) Breuer, M.; Rodi, W. Large eddy simulation of complex turbulent flows of practical interest. In Flow Simulation with HighPerformance Computers II, Notes on Numerical Fluid Mechanics, 2nd Edition; Hirschel, E. H., Ed.; Vieweg Verlag: Braunschweig, Germany, 1996; pp 258−274. (25) 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. (26) Crowe, C. T.; Schwarzkoph, J. D.; Sommerfeld, M.; Tsuji, Y. Multiphase Flows with Droplets and Particles; CRC Press: Boca Raton, FL, 1998. (27) Cundall, P. A.; Strack, O. D. L. A discrete numerical model for granular assemblies. Geotechnique 1979, 9 (1), 47−65. (28) Koch, D. L.; Hill, R. J. Inertial effects in suspension and porousmedia flows. Annu. Rev. Fluid Mech. 2001, 33, 619−647. (29) Issa, R. I. Solution of the implicitly discretised fluid flow equations by operator-splitting. J. Comput. Phys. 1986, 62 (1), 40−65. (30) Li, Y.; Xu, Y.; Thornton, C. A comparison of discrete element simulations and experiments for ‘sandpiles’ composed of spherical particles. Powder Technol. 2005, 160 (3), 219−228. (31) Müller, C. R.; Holland, D. J.; Sederman, A. J.; Scott, S. A.; Dennis, J. S.; Gladden, L. F. Granular temperature: Comparison of Magnetic Resonance measurements with Discrete Element Model simulations. Powder Technol. 2008, 184 (2), 241−253. (32) Müller, C. R.; Scott, S. A.; Holland, D. J.; Clarke, B. C.; Sederman, A. J.; Dennis, J. S.; Gladden, L. F. Validation of a discrete element model using magnetic resonance measurements. Particuology 2009, 7 (4), 297−306. (33) Luo, K.; Yang, S.; Fang, M.; Fan, J.; Cen, K. LES−DEM investigation of the solid transportation mechanism in a 3-D bubbling fluidized bed. Part I: hydrodynamics. Powder Technol. 2014, 256, 385− 394. (34) Yang, S.; Luo, K.; Fang, M.; Fan, J. LES−DEM investigation of the solid transportation mechanism in a 3-D bubbling fluidized bed. 6912

dx.doi.org/10.1021/ie403046q | Ind. Eng. Chem. Res. 2014, 53, 6896−6912