Development of a numerical model for the hydrodynamics simulation of liquid-solid circulating fluidized beds

Powder Technology 348 (2019) 93–104

Contents lists available at ScienceDirect

Powder Technology journal homepage: www.elsevier.com/locate/powtec

Development of a numerical model for the hydrodynamics simulation of liquid-solid circulating fluidized beds Hao Luo a, Chao Zhang a, Jesse Zhu b,⁎ a b

Department of Mechanical and Materials Engineering, The University of Western Ontario, London, ON N6A 3K7, Canada Department of Chemical and Biochemical Engineering, The University of Western Ontario, London, ON N6A 3K7, Canada

a r t i c l e

i n f o

Article history: Received 31 October 2018 Received in revised form 11 February 2019 Accepted 11 March 2019 Available online 12 March 2019 Keywords: Numerical simulation Computational fluid dynamics (CFD) Liquid-solid circulating fluidized bed (LSCFB) Drag model Near wall treatment Multiphase flow

a b s t r a c t Due to extensive applications of the liquid-solid circulating fluidized beds (LSCFBs) in biochemical and petroleum industries, a detailed computational fluid dynamics (CFD) study is crucial to understand the flow characteristics in LSCFBs. In this paper, the hydrodynamics in LSCFBs is numerically investigated using the Eulerian-Eulerian two-phase model combined with the kinetic theory for the granular phase (KTGP). Key factors affecting the simulation results including the drag model, near wall treatment and boundary conditions are investigated, and the CFD model is validated by comparing the numerical results with the experimental data. Among the seven different drag models examined in this study, the adjusted Syamlal O'Brien drag model and the irregular particle drag model are found to provide the best numerical solutions for spherical and irregular particles, respectively. As for the three different near wall treatments tested, the Menter-Lechner near wall treatment is found to provide the best numerical predictions in the near wall region. © 2019 Published by Elsevier B.V.

1. Introduction Liquid-solid fluidization is a process in which a bed of solid particles is suspended in liquid and converted from a solid-like state to a fluidlike state [1]. With unique liquid-solid contacting features, there are numerous advantages to liquid-solid fluidized beds, such as higher contact efficiency and excellent mass and heat transfer. With decades of development, liquid-solid fluidization has been used in different industrial processes, such as biochemical technology, wastewater treatment, petroleum and metallurgical industries [2]. Liquid-solid fluidization can be divided into four regimes. With the increase in the liquid velocity, the fluidization process will go through the fixed bed regime, the conventional fluidization regime, the circulating fluidization regime and the dilute transport regime. For the circulating fluidization regime, all particles can be entrained out of the bed with a high liquid velocity and need to recirculate back to the bottom of the bed. Therefore, a typical LSCFB reactor normally comprises of a riser, a downer, a liquid-solid separator, a top solid-return pipe and a bottom solid-feed pipe [3]. In the recent decades, computational fluid dynamics (CFD) modelling has become an effective tool for investigating the hydrodynamics inside a CFB riser due to the fast development of computer technology and multiphase flow models [4–7]. Generally, there are two main theories for describing multiphase flows: (1) the Eulerian-Lagrangian (E-L) ⁎ Corresponding author. E-mail address: [email protected] (J. Zhu).

https://doi.org/10.1016/j.powtec.2019.03.018 0032-5910/© 2019 Published by Elsevier B.V.

approach where the particulate trajectory model is used for the solidphase and the particle-particle interactions are neglected and (2) the Eulerian-Eulerian (E-E) approach where both phases are treated as continuous phase. In this study, the Eulerian-Eulerian approach is employed since the solid volume fraction in a LSCFB is high and the interactions between particles need to be considered. In the EulerianEulerian approach, the fluid and solid phases are considered as the interpenetrating continuum. The mass and momentum governing equations, which are closed by the constitutive equations, are solved for both phases [8]. The proper CFD modelling of a LSCFB system requires attention on the dynamics of particles flowing in a fluid, turbulent flows and boundary conditions. Among the interactions between particles and fluid, the drag force is the most important one and several drag correlations have been proposed. From the simplest single particle drag correlations [9] to the complex semi-empirical drag correlations, those drag models can be catalogued into three groups: (1) the models based on the pressure drop of the fixed beds, such as the Ergun equation [10] and Gibilaro drag model [11]; (2) the models based on the velocity-voidage correlations, such as the Richardson-Zaki equation [12], Wen-Yu [13], and Syamlal et al. [14] models; and (3) the models based on the energy minimization multi-scale (EMMS) method. In addition, the Gidaspow [15] and Huilin-Gidaspow [16] models were obtained based on the work of Ergun [10] and Wen-Yu [13], and the adjusted Syamlal-O'Brien drag model [17] was proposed to extend the applicability of the SyamalO'Brien model [14] by adjusting the velocity-voidage function parameters to match the experimental minimum fluidization velocity.

94

H. Luo et al. / Powder Technology 348 (2019) 93–104

from Sang [21] where the experiments were conducted using plastic beads with irregular shapes. Despite the numerous studies conducted on multi-phase flows in the literature, none of those studies has comprehensively compared the applicability of different drag models for spherical and irregular particles in the LSCFB systems. And none of the studies have investigated the influence of near wall treatments, boundary conditions and drag models. Therefore, the present work employs the Eulerian-Eulerian CFD model based on KTGP to systematically study the following four aspects: the influence of drag models, near wall treatments for the turbulence model, specularity coefficient, restitution coefficient and the drag model for particles with irregular shapes. The predictions are compared with the experimental data and the agreements are good. 2. Experimental setup of the LSCFB system

Fig. 1. Experimental setup of the LSCFB riser.

For the wall bounded turbulent flows, due to the considerable influence of the walls on the turbulent flows, an accurate representation of the flow in the near wall region should be obtained to ensure accurate numerical solutions. Typically, there are two approaches used for modelling the near-wall region. The first one adopts semi-empirical functions to represent the viscosity affected region instead of solving turbulence equations, such as the standard wall functions and scalable wall functions. The second one modifies the turbulence model to resolve the flow in the near wall region, such as the enhanced wall treatment [18] and Menter-Lechner near wall treatment [18]. For the boundary conditions, the no-slip boundary condition at the wall is suitable for the liquid phase, which is not appropriate for the solid phase. Therefore, Johnson and Jackson [19] proposed a boundary condition, which contains the specularity coefficient and restitution coefficient, to describe the interaction and energy loss between the granular flow and the wall. Furthermore, irregular particles are often used in industrial applications, which significantly affect the hydrodynamic behaviors of the fluidization process. Therefore, studies that quantify the influence of particle shapes in fluidized beds are conducted in this work. By coupling the modified Syamlal-O'Brien drag model [14] and Haider and Levenspiel [20] velocity-voidage function, a new drag model for irregular particles is proposed. The numerical results from the new drag model are validated by comparing them with the experimental data

The experimental data on the liquid-solid two-phase flow in a LSCFB by Razzak [1] and Sang [21] is used in this paper to validate the numerical models. The schematic diagram of the experimental setup is presented in Fig. 1. It consists of two main sections: riser and downer. The riser is made of Plexiglas and is 5.4 m in height and 0.0762 m in diameter. The downer is made of Plexiglas as well and is 5.05 m in height and 0.2 m in diameter. The liquid-solid separator that separates the entrained solids from the liquid is located at the top of the riser and the solids circulation rate measurement device is located near the top of the downer. At the bottom of the riser, there are two liquid distributors, the primary liquid distributors which occupy 19.5% of the cross-sectional area, and the auxiliary liquid distributor which is a porous plate with 4.8% of opening area and controls the recirculating particles flow rate. The particles are injected from the solids feed pipe, by adjusting the auxiliary liquid flow rate, the quantity of recirculating particles from the storage vessel can be controlled. When the auxiliary liquid velocity is zero, there will be no particles entering the riser. Introducing the auxiliary liquid flow to start feeding particles and with the lifting effect from both the auxiliary liquid and primary liquid, all particles can be fluidized and entrained out of the riser. Then they are separated from the liquid in the liquid-solid separator, and ejected into the downer, finally reach to the solids feed pipe again and complete the circulation. In this work, two different operating conditions with spherical glass beads that are used in the experimental work by Razzak [1] are selected to study the effects of different drag models, near wall treatments, specularity coefficient and restitution coefficient, and one operating condition with irregular plastic beads from the experimental work by Sang [21] is chosen to investigate the drag model for irregular particles. Tap water at the ambient temperature is used in all the simulations. The detailed operation conditions and physical properties of the particles and liquid are listed in Table 1. 3. Numerical models The Eulerian-Eulerian based CFD model is used to solve the governing equations for both the continuous liquid phase and discrete solid phase. The governing equation for the continuous liquid phase is closed by the k-ε turbulence model and the governing equation for the discrete solid phase is closed by the models based on kinetic theory of granular

Table 1 Operation conditions and physical properties of the particles and liquid. Parameters

Liquid phase density (kg/m3)

Liquid phase viscosity (kg/(m ⋅ s))

Particle density (kg/m3)

Particle diameter (μm)

Particle sphericity

Superficial liquid velocity (cm/s)

Superficial solid velocity (cm/s)

Operating condition #1 Operating condition #2 Operating condition #3

998.2 998.2 998.2

0.001003 0.001003 0.001003

2500 2500 1520

500 500 580

1 1 0.7

11.2 35 28

0.747 1.193 0.4

H. Luo et al. / Powder Technology 348 (2019) 93–104

3.2.1. Syamlal O'Brien drag model [14] The Syamlal O'Brien drag model gives the correlation for the drag force between the multi-particle system and single particle system as:

phase (KTGP). Those constitutive equations are derived based on different assumptions and empirical correlations. 3.1. Governing equations

C D ð Re; εÞ ¼

The continuity and momentum equations are given as: 
 ∂
* α q ρq þ ∇  α q ρq vq ¼ 0; ∂t

X

αq ¼ 1

95

C Ds ð Re=V r Þ

ð4Þ

Vr2

where Vr is the velocity-voidage function which is originally defined by Richardson-Zaki [12] and modified by Syamlol O'Brien [14] given as:

ð1Þ

q

V r ¼ Re= Res ¼ εnspherical −1


 
 ∂
! ! !2 ! ! α l ρl vl þ ∇  α l ρl vl ¼ −α l ∇p þ ∇  τ l þ α l ρl g þ K sl vs −vl ∂t
 ! !T τ l ¼ α l μ l ∇  vl þ ∇  vl

ð5Þ

The fluid-solid exchange coefficient Kls is defined as: K ls ¼


 
 ∂
! ! !2 ! ! α s ρs vs þ ∇  α s ρs vs ¼ −α s ∇p þ ∇ps þ ∇  τ s þ α s ρs g þ K ls vl −vs ∂t  
 2 ! !T ! τ s ¼ α s μ s ∇  vs þ ∇  v s þ α s λs − μ s ∇  vs I 3

3α s α l ρl 4V 2r ds

 CD

 Re * * vs − vl Vr

ð6Þ

where the drag coefficient is from Dalla Valle [22].   4:8 2 C D ¼ 0:63 þ pffiffiffiffiffiffiffiffiffi ¼ Res

ð2Þ where αq is the volume fraction of phase q. For the continuous liquid phase, a k-ε turbulence model is employed to close the governing equations. Since the dispersed k-ε model is computationally less expensive and predicts the hydrodynamics equally well as the per-phase turbulence model [3], it is used in the simulations and is given as:

Re ¼

4:8 0:63 þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Re=V r

!2

* * ρl ds vs − vl

ð7Þ

μl

The velocity-voidage correlation Vr is obtained from the experiments by Garside and Al-Dibouni [23].

 
  μ t;l ∂ ε  * ðα l ρl ε l Þ þ ∇  α l ρl vl ε l ¼ ∇  α l ∇εl þ α l l C 1ε Gk;l −C 2ε ρl ε l σε kl ∂t

pffiffiffiffiffiffiffipffiffiffiffiffiffiffiffi ∂ ε * ð3Þ ðα l ρl kl Þ þ ∇  α l ρl vl kl −C 2ε l K ls 2kl − 2kl 3Θs kl ∂t  
pffiffiffiffiffiffiffipffiffiffiffiffiffiffiffi μ t;l ∇k þ α l Gk;q −α l ρl εl −K sl 2kl − 2kl 3Θs ¼ ∇  αl σk l

V r −A ¼ 0:06 Res B−V r

ð8Þ

where A ¼ α 4:14 l (

For the solid phase, the KTGP is employed to model the viscosity, stresses and pressure of the solid phase, which are used in the momentum conservation equations. Based on the KTGP, the viscosity, pressure and stresses for the solid phase can be determined by the granular temperature, which is the mean square of a random particle velocity. The constitutive equations for the solid phase are summarized in Table 2.

B¼

c1 α 1:28 l α dl 1

α l ≤ 0:85 c1 ¼ 0:8 α l N 0:85 d1 ¼ 2:65

3.2.2. Gidaspow drag model [15] To cover all flow conditions, Gidaspow [15] combined the Ergun equation [10] and Wen-Yu drag model [13]. When αl N 0.8, the Wen-Yu drag model [13] is adopted:

3.2. Drag models for spherical particles

* * 3 α s α l ρl vs − vl −2:65 αl K sl ¼ C D 4 ds

The drag force is one of the dominate terms in the momentum equation and represents the momentum exchange between phases. Six widely used drag models are investigated in this study.

ð9Þ

Table 2 Constitutive equations for the solid phase. Solids pressure Radial distribution function Solids shear stress Collisional viscosity Kinetic viscosity Frictional viscosity Bulk viscosity Granular conductivity

Ps = αsρsΘs + 2ρs(1 + ess)α2sgO, ssΘs −1 h α s 1=3 i g O;ss ¼ 1−ð Þ α s; max μs = μs, col + μs, kin + μs, fr rffiffiffiffiffiffi 4 Θs μ s;col ¼ α s ρs ds g O;ss ð1 þ ess Þ 5 π pffiffiffiffiffiffiffiffi h i α s ρs ds Θs π 2 μ s;kin ¼ 1 þ ð1 þ ess Þð3ess −1Þα s gO;ss 6ð3 þ ess Þ 5 P s sinϕ μ s;fr ¼ pffiffiffiffiffiffiffi 2 I2D rffiffiffiffiffiffi 4 Θs λs ¼ α 2s ρs ds g O;ss ð1 þ ess Þ 3 π pffiffiffiffiffiffi h i 15ds ρα s Θs π 2 16 1 þ 12 η ð4η−3Þα kΘs ¼ 4ð41−33ηÞ s g O;ss þ 15π ð41−33ηÞηα s g O;ss 5 η¼

Collisional dissipation of energy

γ Θs ¼

12ð1−e2ss Þg O;ss pffiffiffi ρs α 2s Θ3=2 s ds π

1 2 ð1

Lun et al. [22] Ding and Gidaspow [7]

Gidaspow et al. [15] Syamlal et al. [23] Schaeffer [24] Lun et al. [22] Syamlal et al. [23]

þ ess Þ Lun et al. [22]

96

H. Luo et al. / Powder Technology 348 (2019) 93–104

Under the minimum fluidization conditions, Eq. (8) combined with Eq. (5) can be rewritten as

where CD ¼

i 24 h 1 þ 0:15ðα l Res Þ0:687 α l Res

Remf ¼ Remfs

 A þ 0:06B Remfs 1 þ 0:06 Remfs

ð14Þ

When αl ≤ 0.8, the Ergun equation [10] is adopted:

K sl ¼ 150

α s ð1−α l Þμ l 2

α l ds

þ 1:75

* * ρl α s vs − vl ds

ð10Þ

The minimum fluidization condition and terminal settling condition can be considered to be the same for a single particle. Therefore, under the terminal settling condition, we have

C Dts 3.2.3. Huilin-Gidaspow drag model [16] To avoid the discontinuity of the Gidaspow drag model [15] at αl = 0.8, a blending function was used in the Huilin-Gidaspow drag model [16], which is defined as K sl ¼ ψK sl−Ergun þ ð1−ψÞK sl−Wen−Yu

ð11Þ

where ψ¼

1 arctanð262:5ðα s −0:2ÞÞ þ 2 π

K sl ¼

 ρ v* − v* l f s 18 þ 0:33 α s α l−1:8 Re dp

ð12Þ

where the Reynolds number is defined as

Re ¼

ð15Þ

Reynolds number Rets can be calculated by substituting Eq. (7) into Eq. (15). 2
32 pffiffiffiffiffiffiffiffiffiffiffiffiffi0:5 2 −4:87 6 4:8 þ 2:52 4Ar=3 Rets ¼ 4 5 1:26

ð16Þ

where the Archimedes number is related to the drag coefficient and the Reynolds number

3.2.4. Gibilaro drag model [11] Based on the Ergun equation [10], Gibilaro [11] extended its applicability to the dilute particle system by relating the energy dissipation in the fluidized bed with the unrecoverable pressure loss to obtain the particle drag force under the fully expanded condition. The Gibilaro drag model is shown below [11]. 

2 2 3  πdp ρ f U t πdp
¼ ρs −ρ f g 4 2 6

* * α l ρl dp vs − vl

Ar ¼

ρl 3 ds ðρs −ρl Þg μl2

ð17Þ

The minimum fluidization Reynolds number for the multi-particle system Remf can be determined by the experimental data, so, Vr can be obtained. Furthermore, the drag coefficient and velocity-voidage correlation used here are the same as those in the Syamlal-O'Brien drag model [14]. The coefficient A from Eq. (8) is determined based on the experimental data of the bed voidage at the minimum fluidization condition [23], and the relationship between coefficients c1 and d1 is defined as follow to ensure the continuity of B at the condition of α1 = 0.85 [23], where B ¼ c1 α 1:28 ¼ α d11 1

μl Therefore

3.2.5. Single particle drag correlation [9] The liquid-solid fluidization is homogeneous. Hence, the single particle drag correlation is also considered and the results are compared with the multi-particle drag models. The single particle drag correlation used here is from Lewis et al. [9]. 8 24 > > C Ds ¼ > > Re > < 10 C Ds ¼ > > > Re0:5 > > : C Ds ¼ 0:43

Re b 0:4 0:4 b Re b 500

d1 ¼ 1:28 þ

log10 log10

c1

0:85

ð18Þ

For most cases, once the particle properties and the minimum fluidization condition are known, c1 and d1 can be obtained. 3.4. Drag model for irregular particles

ð13Þ

500 b Re b 200000

3.3. Adjustment of the drag model 3.3.1. Adjusted Syamlal-Obrien drag model [17] Since the Syamlal-O'Brien drag model [14] has been widely used, it is often encountered that the operating conditions are remarkably differed from the original experimental conditions. To extend its applicability, the Syamlal-O'Brien drag model [14] can be adjusted by matching the predicted minimum fluidization velocity with the experimental data and the constants c 1 and d 1 can be modified correspondingly.

The drag force between the liquid and solid mainly depends on the local slip velocity, bed voidage, liquid properties, and solid properties, including particle density, size, and shape. However, the six drag models mentioned above are all derived based on spherical particles. It is essential to take the sphericity into consideration in the drag model for irregular particles. Therefore, a new drag model, which is a modified Syamlal-O'Brien drag model, is proposed in this paper. The numerical results from the new drag model are validated by comparing them with the experimental data from Sang [21] where the experiments were conducted using plastic beads with irregular shapes. Since the Syamlal-O'Brien drag model [14] is based on the single particle drag model and the velocity-voidage function for spherical particles, the purpose of the modification is to replace them with the correlations for irregular particles. The single non-spherical particle drag coefficient proposed by Haider and Levenspiel [20] is employed

H. Luo et al. / Powder Technology 348 (2019) 93–104

instead of replacing ε and the turbulent viscosity from separate equations in the near wall region.

in this research, which is given as i 24 h 1 þ f8:1716 expð−4:0655ψÞg  Res ð0:0964þ0:5565ψÞ Res

C Ds ¼

þ

73:69 Res expð−5:0748ψÞ Res þ 5:378 expð6:2122ψÞ

ð19Þ

The shape factor ψ ¼ where s is the surface area of a sphere which has the same volume as the irregular particle, and S is the actual surface area of the irregular particle [24]. The original velocity-voidage function for spherical particles from Richardson-Zaki [12] is modified by Cleasby and Fan [25] in which a function for the shape factor is included in “n” value to account for the effects of irregular shapes. s S,

V r ¼ εn−1

ð20Þ

where n = nsperical × (ψ)γ, γ = − 2.9237 Re t−0.363(ψ)0.884 nspherical is the “n” value form RZ equation [12] given as 8 nsperical > > >n > < sperical nsperical > > > > : nsperical nsperical

¼ 4:65 þ 20d=D ¼ ð4:4 þ 18d=DÞ Ret −0:03 ¼ ð4:4 þ 18d=dÞ Ret ¼ 4:4 Ret −0:1 ¼ 2:4

Ret ≤ 0:2 0:2 b Ret ≤ 1 1 b Ret ≤ 200 200 b Ret ≤ 500 Ret N 500

Therefore, substituting the correlations for irregular particles, Eqs. (19) and (20), into Eq. (4), the drag coefficient for irregular particles can be obtained as the following. CD ¼

(  0:0964þ0:5565ψ ) 24 Re 1 þ ½8:1716 expð−4:0655ψÞ  ReVr Vr þ

73:69 Re expð−5:0748ψÞ

97

ð21Þ

Vr2 Re þ 5:378Vr3 expð6:2122ψÞ

where Vr is defined by Eq. (20).

4. Numerical methodology To simulate the two-phase flows in the LSCFB shown in Fig. 1, the riser is simplified to a 2D-planar as shown in Fig. 2 and the mesh information can be seen from Table 3. Due to the potential instantaneous non-axisymmetric flow structures within LSCFBs, the 2D planar mesh is created and the results are time averaged for 20s after reaching a stable condition. In addition, with the use of the enhanced wall treatment, Y+ ≤ 1 should be satisfied for the first grid adjacent to the wall. Therefore, the finer grid is used in the near wall region with an expansion ratio of 1.05 for the cell size from the wall to the center of the bed. Additionally, to correctly represent the complex flow structures at the inlet, the mesh in the inlet region has been refined as shown in Fig. 3 and the expansion ratio is set as 1.05 as well. The mesh independence is examined for operating condition #3 using three different grids, 60 × 1500, 100 × 2500 and 120 × 2500, in the x and y directions, respectively. The radial profiles of the solids concentration at different bed heights, H = 1.01 m, 2.02 m, 3.03 m and 3.82 m, are compared. It is found that the difference in the results for all three meshes is b0.5%. Therefore, the medium mesh (100 × 2500) is used to reduce the calculation time while ensuring accurate results. The mesh information for different operating conditions is given in Table 3. For the boundary conditions, at the inlet, which is located at the bottom of the riser, both the liquid and particles are of uniform velocities. At the outlet, outflow condition is used due to the fully developed flow condition at the outlet. On the wall, the no-slip condition is used for the liquid phase and the partial slip Johnson and Jackson [19] boundary condition is used for the solid phase. The dispersed k-ε turbulence model is used for the liquid phase while the particle-particle collision restitution coefficient is set as 0.95 for the solid phase. The phase coupled SIMPLE scheme is used for the pressure-velocity coupling.

3.5. Near wall treatment In this paper, three different near wall treatments including the scalable wall function [18], the enhanced wall treatment [18] and the Menter-Lechner near wall treatment [18] are implemented and examined. 3.5.1. Scalable wall functions [18] The scalable wall function is modified based on the standard wall function to avoid the computational deterioration when y* b 11. It adds a selector as shown below. ye ¼ MAX ðy ; ylimit Þ

ð22Þ

where y∗limit =11.225. Hence, if y* N 11.225, the standard and scalable wall functions are identical. It should be noted that the y* and Y+ are approximately equal in the equilibrium turbulent boundary layer. 3.5.2. Enhanced wall treatment [18] By modifying the turbulence model, the enhanced wall treatment ensures to resolve the viscous sublayer with a refined mesh where the first wall-adjunct grid meets Y+ = 1. The two-layer approach is employed to specify both ε and the turbulent viscosity in the nearwall cells [18]. 3.5.3. Menter-Lechner ε-equation [18] The Menter-Lechner near wall treatment is introduced to avoid the drawbacks of the enhanced wall treatment. In the Menter-Lechner wall treatment, a source term is added to the k − ε turbulence model

Fig. 2. Schematics of the LSCFB riser.

98

H. Luo et al. / Powder Technology 348 (2019) 93–104

Table 3 Mesh information for different operating conditions. Parameters

Domain size (m)

Number of control volumes

Wall space for first grid (m)

Expansion ratio along the radius

Expansion ratio along the axis

Maximum aspect ratio

Operating condition #1 Operating condition #2 Operating condition #3

0.0762 × 5.97 0.0381 × 5.97 0.0762 × 5.2

100 × 2500 100 × 2500 100 × 2500

0.00015 0.000052 0.000064

1.05 1.05 1.05

1.05 1.05 1.05

16.73 41.12 29.48

The power law is chosen to discretize the convection terms for the k-ε turbulence model and granular temperature while the QUICK is chosen for the mass and momentum governing equations. Furthermore, the time step size is set as 1 × 10−04 s and the convergence criteria is set as 5 × 10−05. The parameters in the proposed drag model for irregular particles can be determined based on the sphericity of the particles used in the experimental work by Sang [21], which is 0.7. The proposed drag model is compiled into Fluent solver by the User Defined Function (UDF).

5. Results and discussion The numerical models described above are employed to predict the flow field and hydrodynamics of a LSCFB riser. The effects of the drag models for spherical and irregular particles, the near wall treatments for the turbulence flow and the coefficients for the Johnson and Jackson [19] boundary conditions are examined by comparing the numerical results with the experimental data. In addition, the effects of some critical parameters affecting the predictions are analyzed.

Fig. 3. Computational domain and mesh.

5.1. Drag models for spherical particles To investigate the influence of drag models on the numerical results for risers using spherical particles, the simulations are conducted using six different drag correlations with the enhanced wall treatment under two operating conditions, Operating Conditions #1 and #2 as shown in Table 1. For the adjusted Syamlal-O'Brien drag model [17], the parameters for the velocity-voidage correlation c1 and d1, which only depend on the properties of particles and fluid, are adjusted to 0.304 and 8.605, respectively, for both operating conditions to match with the experimental data. The comparisons for the radial solids holdup profiles using different drag modes with the enhanced wall treatment are shown in Figs. 4 and 5 under two operating conditions, respectively. All the predicted radial solids holdup distributions have the same trend. That is due to the frictions and no slip conditions between the liquid and the wall. The liquid velocity decreases towards the wall and finally approaches 0 at the wall,

Fig. 4. Comparisons of the radial distributions of the solids holdup using different drag models under Ul = 11.2 cm/s and Us = 0.747 cm/s. (a) H = 1.01 m and (b) H = 3.82 m.

H. Luo et al. / Powder Technology 348 (2019) 93–104

Fig. 5. Comparisons of the radial distributions of the solids holdup using different drag models under Ul = 35 cm/s and Us = 1.193 cm/s. (a) H = 1.01 m and (b) H = 3.82 m.

which causes to the same trend for the solid velocity. Therefore, the solids concentration is higher at the near wall region and lower at the central region. However, all the simulations tend to overestimate the solids holdup compared to the experimental data and the increase in the solids holdup towards the wall is not as obvious as that shown in the experiment. The difference between the numerical results and experimental data are listed in Tables 4 and 5 under different operating conditions. It is clear from Tables 4 and 5, for the central region, among different drag models, the adjusted Syamlal O'Brien drag model [17] gives the best result while the Gibilaro model [11] results in the worst result. The results from the Gidaspow [15] and Huilin-Gidaspow [16] drag models are almost the same, which are better than that from the Gibilaro model [11], but worse than that from the Syamlal O'Brien drag model [14] as shown in Tables 4 and 5. It is noticed that the

99

Gidaspow [15] and Huilin-Gidaspow [16] drag models give smaller errors in the near wall region than other models. However, those two models give large errors in the center region as shown in Tables 4 and 5. In summary, the adjusted Syamlal O'Brien drag model [17] gives the best overall results. However, the increasing trend for the solids holdup towards to the wall is not correctly predicted by all those models, i.e., they predicted more flat radial profiles for the solids holdup than those from the experimental data. The granular temperature is also an important parameter. By applying the kinetic theory of the granular phase, the solids pressure and solids viscosity can be determined by the granular temperature. Generally, the granular temperature represents the random kinetic energy of particles per unit mass and it greatly depends on the particle velocity fluctuations. Therefore, a higher granular temperature means a higher velocity fluctuation and more intense particle collisions, which results in a lower solids velocity and higher solids concentration. Figs. 6 and 7 show the radial distribution of the granular temperature under two different operating conditions. It is obvious that the predicted granular temperature from the adjusted Syamlal O'Brien drag model [17] is lower than that from the Syamlal O'Brien drag model [14]. The granular temperatures from the Gidaspow [15] and Huilin-Gidaspow [16] drag models are almost identical and higher than that from the Syamlal O'Brien model [14]. Besides, the Gibilaro drag model [11] provides the highest radial granular temperature distribution and results in the highest solids concentration. From the comparisons for both solids holdups and granular temperatures using different drag models, it can be seen among all the drag models, the Gibilaro model [11] tends to provide the least accurate results. It is probably due to the fact that it was developed from the Ergun equation [10], which was based on the pressure drops of different fixed beds. Even it was modified by considering the fully expanded bed condition, it still cannot correctly represent the flow filed in LSCFBs. The second least accurate drag models are Gidaspow model [15] and HuilinGidaspow model [16]. Those two models give almost identical results and were developed based on the semi-empirical correlation of Ergun equation [10] and Wen-Yu model [13]. For the LSCFB under the two operating conditions considered in this study, the bed voidage is always N0.85, hence, the Wen-Yu drag model [13] is used in the Gidaspow and Huilin-Gidaspow models [15,16]. Therefore, it is found the WenYu drag model [13] performs better than Gibilaro drag model [11], but still needs to be improved to be applied in LSCFBs. As for the Syamlal O'Brien drag model [14], the derivation of the semi-empirical correlation was based on both gas-solid and liquid-solid fluidized beds, therefore, it can provide better predictions for the LSCFBs. However, the general drag models cannot always provide satisfactory predictions for different cases. It is shown in Figs. 4 and 5, by adjusting the parameters c1 and d1 correspondingly based on the case-specified minimum fluidization conditions, the adjusted Syamlal-O'Brien drag model [17] can provide better predictions. In addition, the single particle correlation [9] can provide relatively satisfactory results, it is mainly due to the low solids concentration and relative homogeneous flow structures in the LSCFB, which is similar to the dilute single particulate flow. However, beside the improvement of the solids holdup distribution along the radial direction, it is noticed that the adjusted Syamlal O'Brien

Table 4 Difference between the numerical results and experimental data using different drag models at H = 3.82 m under Ul = 11.2 cm/s and Us = 0.747 cm/s. r/R

Experimental solids holdup

0 0.0895 0.2 0.0895 0.49 0.0916 0.64 0.0936 0.76 0.1017 0.86 0.1096 0.95 0.1138 Averaged difference

Gibilaro

Syamlal O'Brien

Adjusted Syamlal O'Brien

Gidaspow

Huilin-Gidaspow

56.42% 56.42% 53.38% 50.64% 39.13% 30.02% 28.30% 44.90%

10.85% 10.85% 8.31% 5.99% 2.45% 9.48% 9.49% 8.20%

2.23% 2.23% 3.93% 5.98% 11.50% 17.88% 17.40% 8.74%

22.73% 22.75% 20.02% 17.55% 8.39% 1.09% 0.31% 13.26%

24.49% 24.50% 21.48% 18.74% 9.12% 1.40% 1.05% 14.40%

100

H. Luo et al. / Powder Technology 348 (2019) 93–104

Table 5 Difference between the numerical results and experimental data using different drag models at H = 1.01 m under Ul = 35 cm/s and Us = 1.193 cm/s. r/R

Experimental solids holdup

0 0.036 0.2 0.036 0.49 0.037 0.64 0.0379 0.76 0.0395 0.86 0.0437 0.95 0.0471 Averaged difference

Gibilaro

Syamlal O'Brien

Adjusted Syamlal O'Brien

Single particle correlation

Huilin-Gidaspow

22.03% 22.08% 20.46% 20.18% 18.03% 10.59% 5.99% 17.05%

8.83% 8.92% 7.51% 6.83% 4.46% 3.02% 5.52% 6.44%

5.56% 6.39% 4.86% 3.96% 1.34% 6.18% 8.70% 5.28%

10.92% 11.03% 9.76% 8.89% 6.46% 1.14% 4.06% 7.47%

11.33% 11.39% 9.89% 9.47% 7.24% 0.25% 2.91% 7.50%

model [17] works well for all the conditions except for the one at the height = 3.82 m as shown in Figs. 4 and 5 because the solids suspension is too dilute in the upper region of the riser. Though, it is the one that gets the closest to the experimental data than the other drag models, which also highlights the need of further improvement of the drag model that can predict the accurate solid holdup distribution along the axial direction. 5.2. Near wall treatments

the improvement in the numerical predictions. Furthermore, from the comparison of the granular temperature shown in Fig. 9, there is a modest increase in the granular temperature near the wall by using the Menter-Lechner near wall treatment when compared with that from the enhanced wall treatment. This results in a higher velocity fluctuation and more intense particle collisions, and leads to a lower solids velocity and higher solids concentration near the wall. As for the scalable wall function, it is clear the distributions of both solids holdup and granular temperature are almost flat, which does not agree with the experimental data. Therefore, the scalable wall function is not suitable for the fine mesh (Y+ is b1) used in this study. To evaluate the influence of the near wall treatment on the performance of different drag models, the simulations using different drag models with the Menter-Lechner near wall treatment [18] are carried out in this study. The comparison of the radial solids holdup is shown in Fig. 10 and the difference between numerical results and experimental data is shown in Table 7. It is noticed that the trend of the performance of different drag models using the Menter-Lechner near wall treatment [18] is the same as that shown in Fig. 5 where the enhanced wall treatment was applied, i.e. the adjusted Syamlal-O'Brien drag model still provides the best numerical predictions compared with other drag models. However, it is observed the solids holdup at the near wall region is significantly improved by using the MenterLechner near wall treatment. In addition, it is clear as shown in Table 7 that using the adjusted Syamlal O'Brien drag model and Menter-Lechner near wall treatment will result in the best agreement between the numerical prediction and experimental data.

Based on the comparisons between different drag models, it can be seen the predicted solids holdup at the central region can be improved by using the adjusted Syamlal O'Brien drag model [17]. However, this model gives higher errors at the near wall region, which leads to investigations on the near wall treatments. For a wall bounded flow, an accurate representation at near wall region can significantly improve the numerical results. Therefore, three different near-wall treatment methods, the scalable wall function [18], enhanced wall treatment [18] and Menter-Lechner method [18], are investigated in this study. The comparisons for the radial profiles of the solids holdup and granular temperature are illustrated in Figs. 8 and 9 under Operating Condition #2, and the difference between the numerical results and experimental data is presented in Table 6. It can be seen from Fig. 8 and Table 6 that the Menter-Lechner near wall treatment gives the best agreement with the experimental data at the near wall region. The scalable near wall function results in an almost flat solids distribution, i.e. no increase in the solids holdup at the near wall region. With the adjusted Syamlal O'Brien drag model [17], the solids holdup distributions at the central region using all near wall treatments show good agreements with the experimental data, while the difference between numerical results and experimental data at the near wall region is significantly smaller by using the Menter-Lechner near wall treatment than other near wall treatments, which represents

The boundary conditions are very important for the wall bounded flows, such as the multiphase flows in LSCFBs. For the liquid phase, the near wall treatment is used to resolve the flow in the near wall

Fig. 6. Comparison of the radial distributions of the granular temperature using different drag models at H = 3.82 m under Ul = 11.2 cm/s and Us = 0.747 cm/s.

Fig. 7. Comparison of the radial distributions of the granular temperature using different drag models at H = 1.01 m under Ul = 35 cm/s and Us = 1.193 cm/s.

5.3. Specularity and restitution coefficients

H. Luo et al. / Powder Technology 348 (2019) 93–104

101

Table 6 Difference between the numerical results and experimental data using different near wall treatments at H = 1.01 m under Ul = 35 cm/s and Us = 1.193 cm/s. r/R

Experimental solids holdup

0 0.036 0.2 0.036 0.49 0.037 0.64 0.0379 0.76 0.0395 0.86 0.0437 0.95 0.0471 Averaged difference

Fig. 8. Comparison of the radial distributions of the solids holdup using different near wall treatments at H = 1.01 m under Ul = 35 cm/s and Us = 1.193 cm/s.

region. While for the solid phase, the interactions between the particles and the wall are modeled using the specularity coefficient and restitution coefficient, which are discussed in this section. The restitution coefficient stands for the ratio of the velocity change during the collision between the particles and the wall, which is implemented in the Johnson-Jackson granular boundary conditions [19]. It varies from 0 to 1, from the inelastic collision to the elastic collision. It was found the restitution coefficient is close to unity for particles and the wall in LSCFBs due to the lubrication effect brought by liquid film [26]. Therefore, the restitution coefficients of 0.85, 0.9, 0.95 and 0.99 are tested in this study. Similarly, the specularity coefficient is also specified in the Johnson-Jackson boundary conditions [19]. When it is zero, the condition is described as zero shear at the wall, while the unity represents that there is a significant amount of lateral momentum transfer at the wall [18]. To investigate the sensitivity of predicted hydrodynamics to the specularity coefficient, the specularity coefficient of 0.00005, 0.0005, 0.005, 0.05 and 0.5 are selected. The comparisons of radial solids holdup distributions using different restitution and specularity coefficients are shown in Figs. 11 and 12. It can be seen from Fig. 11 that there is no notable variation for the radial solids holdup distributions under different restitution coefficients at both locations along the riser, which indicates the flow field in the LSCFB is not sensitive to the restitution coefficient between the particles

Fig. 9. Comparison of the radial distributions of the solids granular temperature using different near wall treatments at H = 1.01 m under Ul = 35 cm/s and Us = 1.193 cm/s.

Enhanced wall treatment

Menter-Lechner near wall treatment

Scalable wall function

6.22% 6.31% 4.95% 3.93% 1.34% 6.20% 8.73% 5.38%

5.56% 5.56% 4.35% 4.85% 4.00% 1.60% 2.34% 4.04%

6.17% 6.19% 3.30% 0.92% 3.16% 12.49% 18.87% 7.30%

and wall. Fig. 12 shows the comparison of the radial solids holdup under different specularity coefficients. It can be clearly seen that there is not much of a difference when the specularity coefficient is between 0.00005 and 0.05. However, there is a noticeable increase in the solids holdup at the near wall region when the specularity coefficient is 0.5, which means within the LSCFB system, the solids distribution is not sensitive to the specularity coefficient until it reaches a critical value. However, the specularity coefficient at 0.5 is not physically possible since the particle-wall collisions in LSCFBs should be close to elastic. Therefore, the specularity coefficient and restitution coefficient are chosen as 0.0005 and 0.95 for the rest of simulations. At last, the numerical solutions with the adjusted Syamlal-O'Brien drag model [17] and the Menter-Lechner near wall treatment under two different operating conditions are shown in Figs. 13 and 14. It can be seen the simulation results under both operating conditions at different heights in the riser have good agreements with the experimental data, i.e. lower solids concentration at the central region and higher solids concentration near the wall. However, with the increase of the height along the riser, the agreement between the numerical results and experimental data is not as strong as those at the lower part of the riser. 5.4. Drag models for irregular particles Since irregular particles are often used in industrial applications and the sphericity of the particles will affect the hydrodynamics within the LSCFB system, it is essential to take the sphericity into the consideration during the design or scale-up of a fluidized bed. In the present work, to investigate the influence of the particle shapes on the drag force, the

Fig. 10. Comparison of the radial distributions of the solids holdup using different drag models at H = 1.01 m under Ul = 35 cm/s and Us = 1.193 cm/s.

102

H. Luo et al. / Powder Technology 348 (2019) 93–104

Table 7 Difference between the numerical results and experimental data using different drag models at H = 1.01 m under Ul = 35 cm/s and Us = 1.193 cm/s. r/R

Experimental solids holdup

0 0.036 0.2 0.036 0.49 0.037 0.64 0.0379 0.76 0.0395 0.86 0.0437 0.95 0.0471 Averaged difference

Adjusted Syamlal-O'Brien

Syamlal-O'Brien

Gidaspow

Gibilaro

5.56% 5.56% 4.35% 4.85% 4.00% 1.60% 2.34% 4.04%

8.67% 8.83% 6.41% 6.78% 6.10% 0.64% 1.04% 5.50%

11.36% 11.36% 9.32% 9.82% 9.04% 3.46% 0.81% 7.88%

22.00% 22.00% 19.57% 19.82% 19.27% 14.03% 8.81% 17.93%

simulations are carried out by employing the modified drag model for irregular particles given by Eq. (20) with the Menter-Lechner near wall treatment and Johnson-Jackson boundary conditions under the Operating Condition #3 in Table 1. The numerical results are compared with the experimental data. In addition, the constants in the adjusted Syamlal O'Brien model [17] are adjusted for irregular particles, i.e. c1 and d1 are adjusted to 0.282 and 9.074, respectively, for Operating Condition #3. The comparison for radial solids holdup profiles using different drag models is shown in Fig. 15. It is clear the radial solids holdup distribution for irregular plastic beads has the analogous manner as spherical glass beads, i.e. lower solids concentration at the central region and higher solids concentration near the wall, which indicates different particle shapes will not alter the basic hydrodynamic behaviors of the LSCFB system. It can be seen that all the simulations tend to overestimate the solids holdup compared to the experimental data. Among those different drag models, it is evident the irregular particle drag model gives the best results, the adjusted Syamlal O'Brien model [17] is the next best one, and followed by the Gidaspow model [15]. Additionally, due

Fig. 11. Comparisons of the radial distributions of the solids holdup using different restitution coefficients for the specularity coefficient of 0.0005 under Ul = 35 cm/s and Us = 1.193 cm/s.

Fig. 12. Comparisons of the radial distributions of the solids holdup using different specularity coefficient for the restitution coefficient of 0.95 under Ul = 35 cm/s and Us = 1.193 cm/s.

Fig. 13. Comparisons of the radial distributions of the solids holdup under Ul = 11.2 cm/s and Us = 0.747 cm/s at different axial locations.

H. Luo et al. / Powder Technology 348 (2019) 93–104

103

Fig. 16. Comparison of the radial distributions of the granular temperature for irregular particles using different drag models at H = 3.98 m under Ul = 28 cm/s and Us = 0.4 cm/s.

Fig. 14. Comparisons of the radial distributions of the solids holdup under Ul = 35 cm/s and Us = 1.193 cm/s at different axial locations.

to the use of the Menter-Lechner near wall treatment, the radial solids concentration near the wall shows good agreement with the experimental data. The radial granular temperature distributions using different drag models under Operating Condition #3 are shown in Fig. 16. It is evident that the granular temperature predicted by the irregular particle drag model is lower than the adjusted Syamlal O'Brien model [17] drag model and the Gidaspow [15] drag model gives the highest granular temperature distribution. The predictions of the granular temperature are consistent with the solids holdup distributions. It can be seen among all the drag models, the Gidaspow drag model [15] always gives the least accurate results. It is reasonable as mentioned in the previous section, when the bed voidage is N0.85, the Wen-Yu drag model [13] is used and it does not take the shape factor

into consideration. Therefore, the Gidaspow drag model [15] is not suitable for the irregular LSCFB system. As for the adjusted Syamlal O'Brien drag model [17] and the irregular particle drag model, they are obtained by the same concept, which is by including the velocity-voidage function in the single particle drag correlation to obtain the drag model for multi-particle systems. For the adjusted Syamlal O'Brien drag model [17], both the single particle drag model from Dalla Valle [22] and velocity-voidage function from Garside and Al-Dibouni [23] were obtained based on spherical particles in gas or liquid-solid systems. Therefore, by modifying c1 and d1 which both depend on particle properties, there is an improvement for the adjusted Syamlal O'Brien drag model [17] compared with the Gidaspow drag model [15]. However, for the irregular particle drag model, the sphericity has been directly taken into consideration in both single particle drag model and velocity-voidage function. The single non-spherical particle drag model is adopted from Haider and Levenspiel [20] which was obtained from 419 isometric data points for Re b 2.5 × 105 (including irregular particles like octahedrons, cubes, tetrahedrons, disks, etc.) and 408 data points for Re b 2.6 × 105 (including spherical particles to test the applicability when the sphericity is 1). Then, the velocity-voidage function is chosen from Cleasby and Fan [25], which was obtained from the experimental studies for irregular particles such as sand, anthracite and flints, etc. Therefore, by directly taking sphericity into consideration, the irregular particle based semi-empirical drag model is more accurate than the adjusted Syamlal O'Brien drag model [17], which was based on the spherical particles.

6. Conclusions

Fig. 15. Comparison of the radial distributions of the solids holdup for irregular particles using different drag models at H = 3.98 m under Ul = 28 cm/s and Us = 0.4 cm/s.

A numerical study on the effects of the drag models, near-wall treatments and wall boundary conditions for the predictions of the turbulent liquid-solid two-phase flows in a fluidized bed has been carried out. For the spherical particle systems, it is found that the adjusted SyamlalO'Brien drag model [17] provides the best agreement with the experimental data at the central region while the Menter-Lechner near wall treatment [18] gives a more realistic solution at the near wall region. In addition, the numerical predictions are not sensitive to the specularity coefficient and restitution coefficient for the liquid-solid two-phase flows in a fluidized bed. For the irregular particle systems, it is concluded that by including the shape factor into the single non-spherical particle drag model and velocity-voidage function to account for the effect of the irregular particles, the prediction accuracy can be improved. In future works, the comprehensive numerical model proposed in this study will be adopted to investigate the influence of the superficial

104

H. Luo et al. / Powder Technology 348 (2019) 93–104

liquid velocity, superficial solid velocity and particle density acts on the hydrodynamics of LSCFBs. References [1] S.A. Razzak, S. Barghi, J.-X. Zhu, Axial hydrodynamic studies in a gas–liquid–solid circulating fluidized bed riser, Powder Technol. 199 (2010) 77–86, https://doi.org/10. 1016/J.POWTEC.2009.05.014. [2] A. Atta, S.A. Razzak, K.D.P. Nigam, J.-X. Zhu, (Gas)−liquid−solid circulating fluidized bed reactors: characteristics and applications, Ind. Eng. Chem. Res. 48 (2009) 7876–7892, https://doi.org/10.1021/ie900163t. [3] A. Dadashi, J. Zhu, C. Zhang, A computational fluid dynamics study on the flow field in a liquid–solid circulating fluidized bed riser, Powder Technol. 260 (2014) 52–58. [4] W. Wang, B. Lu, N. Zhang, Z. Shi, J. Li, A review of multiscale CFD for gas–solid CFB modeling, Int. J. Multiph. Flow. 36 (2010) 109–118, https://doi.org/10.1016/J. IJMULTIPHASEFLOW.2009.01.008. [5] L. Ratschow, R. Wischnewski, CFD-simulation of a circulating fluidized bed riser, Particuology 7 (2009) 283–296, https://doi.org/10.1016/J.PARTIC.2009.04.005. [6] J.L. Sinclair, R. Jackson, Gas-particle flow in a vertical pipe with particle-particle interactions, AICHE J. 35 (1989) 1473–1486, https://doi.org/10.1002/aic.690350908. [7] J. Ding, D. Gidaspow, A bubbling fluidization model using kinetic theory of granular flow, AICHE J. 36 (1990) 523–538, https://doi.org/10.1002/aic.690360404. [8] H. Enwald, E. Peirano, A.-E. Almstedt, Eulerian two-phase flow theory applied to fluidization, Int. J. Multiph. Flow. 22 (1996) 21–66, https://doi.org/10.1016/S03019322(96)90004-X. [9] D. Kunii, O. Levenspiel, D. Kunii, O. Levenspiel, Entrainment and elutriation from fluidized beds, J. Chem. Eng. Jpn. 2 (1969) 84–88, https://doi.org/10.1252/jcej.2.84. [10] M. Annaland, A. Bokkers, M. Goldschmidt, O.O. Olaofe, M.A. van der Hoef, H. Kuipers, Development of a Multi-Fluid Model for Poly-Disperse Dense Gas–Solid Fluidised Beds, Part II: Segregation in Binary Particle Mixtures, 2009, https://doi.org/10. 1016/j.ces.2009.06.043. [11] L.G. Gibilaro, R. Di Felice, S.P. Waldram, P.U. Foscolo, Generalized friction factor and drag coefficient correlations for fluid-particle interactions, Chem. Eng. Sci. 40 (1985) 1817–1823, https://doi.org/10.1016/0009-2509(85)80116-0.

[12] J.F. Richardson, W.N. Zaki, Sedimentation and fluidisation: part I, Chem. Eng. Res. Des. 75 (1997), https://doi.org/10.1016/S0263-8762(97)80006-8 S82–S100. [13] Y.H. Wen, C.Y. Yu, Mechanics of fluidization, Chem. Enginnering Prgress Symp. Ser 1966, pp. 100–111. [14] M. Syamlal, T. O'Brien, The Derivation of a Drag Coefficient Formula from VelocityVoidage Correlations, 1994. [15] D. Gidaspow, in: D.B.T.-M. F, F. Gidaspow (Eds.), 9 - Kinetic Theory Approach, Academic Press, San Diego 1994, pp. 239–296, https://doi.org/10.1016/B978-0-08051226-6.50013-3. [16] L. Huilin, D. Gidaspow, Hydrodynamics of binary fluidization in a riser: CFD simulation using two granular temperatures, Chem. Eng. Sci. 58 (2003) 3777–3792, https://doi.org/10.1016/S0009-2509(03)00238-0. [17] M. Syamlal, T.J. O'Brien, Fluid dynamic simulation of O3 decomposition in a bubbling fluidized bed, AICHE J. 49 (2003) 2793–2801, https://doi.org/10.1002/aic. 690491112. [18] Ansys Fluent 15.0, theory guide, November 2013. [19] P.C. Johnson, R. Jackson, Frictional–collisional constitutive relations for granular materials, with application to plane shearing, J. Fluid Mech. 176 (1987) 67, https://doi. org/10.1017/S0022112087000570. [20] A. Haider, O. Levenspiel, Drag coefficient and terminal velocity of spherical and nonspherical particles, Powder Technol. 58 (1989) 63–70, https://doi.org/10.1016/ 0032-5910(89)80008-7. [21] L. Sang, Particle Fluidization in Upward and Inverse Liquid-Solid Circulating Fluidized Bed, Univ. West. Ontario, London, 2013. [22] J.M. DallaValle, A. Klemin, Micromeritics: The Technology of Fine Particles, vol. xiv, 1943 , (428 p) file://catalog.hathitrust.org/Record/001512807. [23] J. Garside, M.R. Al-Dibouni, Velocity-voidage relationships for fluidization and sedimentation in solid-liquid systems, Ind. Eng. Chem. Process. Des. Dev. 16 (1977) 206–214, https://doi.org/10.1021/i260062a008. [24] H. Wadell, Sphericity and roundness of rock particles, J. Geol. 41 (1933) 310–331, http://www.jstor.org/stable/30058841. [25] J.L. Cleasby, K.S. Fan, Fluidization and Expansion of Filter Media, 1981. [26] Y. Cheng, J. Zhu, CFD modelling and simulation of hydrodynamics in liquid-solid circulating fluidized beds, Chem. Eng. Sci. 63 (2008) 3201–3211.