Comparative analysis of heterogeneous gas-solid flow using dynamic cluster structure-dependent drag model in risers

International Journal of Multiphase Flow 122 (2020) 103126

Contents lists available at ScienceDirect

International Journal of Multiphase Flow journal homepage: www.elsevier.com/locate/ijmulflow

Comparative analysis of heterogeneous gas-solid flow using dynamic cluster structure-dependent drag model in risers Jiang Xiaoxue a, Li Dan a, Wang Shuyan b,∗, M. Hassan a, Cai Wenjian a, Chen Weiqi a, Lu Huilin a,∗ a b

School of Energy Science and Engineering, Harbin Institute of Technology, Harbin, 150001, China School of Petroleum Engineering, Northeast Petroleum University, Daqing, 163318, China

a r t i c l e

i n f o

Article history: Received 20 June 2019 Revised 7 July 2019 Accepted 27 September 2019 Available online 28 September 2019 Keywords: Dynamic cluster structure-dependent drag model Two-fluid model Kinetic theory of granular flow Minimization of energy dissipation rate Heterogeneous structure Fluidized bed

a b s t r a c t An energy dissipation minimization-based dynamic cluster structure-dependent (DCSD) drag model is proposed using Euler-Euler two-fluid model (TFM) and kinetic theory of granular flow (KTGF). The DCSD drag model consists of a set of transient nonlinear equations which includes eight balance equations and one extreme value equation of the minimization of energy dissipation rate. The criterion of cluster existence is proposed to identify the heterogeneous flow and homogenous flow of gasparticles suspension. The intermittency factor is defined to describe the occurrence time of clusters. Three clustering mechanisms named collision-dominant (CD), collision-hydrodynamic-dominant (CHD) and hydrodynamic-dominant (HD) are identified as a result of the compromise in competition of energy dissipation rate components of hydrodynamic interactions and interactions of collision of particles. The CHD for particle clustering is the central mechanism, and CD is next in importance to HD. The predicted cluster diameter and solids volume fraction of clusters are compared to calculations using empirical correlations in literature. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Circulating fluidized bed (CFB) risers are widely used in coal combustion and gasification, oil refining and several other industrial applications. Experiments and numerical simulations showed that the existence of particle clusters effects mass and heat transfer and chemical reactions in CFB risers (Yerushalmi et al., 1976; Grace and Tuot, 1979; Molerus et al., 1997; Basu, 2015; Yates et al., 2016). The identification of particles clusters has been studied using statistical methodology from experimental measurements and simulation results (Grace et al., 1997; Gidaspow et al., 2009; Arastoopour et al., 2016). The findings of measured clusters have been reviewed by Cahyadi et al. (2017) in CFB risers. With the progression of cluster studies, it is found that particle clustering causes significant variation of momentum transfer between the gas and solid phases, making the numerical simulation of flow of gas and particles extremely difficult by means of Euler-Euler two-fluid model (TFM) with the kinetic theory of granular flow (KTGF) (Zhang et al., 2002; Breault et al., 2006). Syamlal & O’Brien (2003) proposed a drag model to explain the

∗

Corresponding authors. E-mail addresses: [email protected] (W. Shuyan), [email protected] (L. Huilin). https://doi.org/10.1016/j.ijmultiphaseflow.2019.103126 0301-9322/© 2019 Elsevier Ltd. All rights reserved.

higher than expected drag coefficients measured in CFBs. A scaled drag model using TFM is implemented into simulations of fluid catalytic cracking (FCC) particles as a function of scale factor and solid volume fraction according to the Wen & Yu drag correlation (McKeen and Pugsley, 2003; Li et al., 2008). The value of the scale factor is in the range of 0.2–0.3. These empirical drag models are relatively simple, and are easy to implement in computational fluid dynamics (CFD) codes. There are at least two approaches, the filtered multiscale models and the heterogeneous structure-based drag models, used for modeling flow of heterogeneous structures of clusters and bubbles in fluidized beds. The filtered multiscale models include the filtered two-fluid model (fTFM) (Agrawal et al., 2001), filtered drag model (Parmentier et al., 2012), filtered drag-stress model (Scherabueaer et al., 2014), anisotropic filtered model (Cloete et al., 2018), enhanced filtered model (Gao et al., 2018) and filtered stress model (Mouallem et al., 2019). The conservation equations of the fTFM are derived by applying space filtering operators, assuming the flow of heterogeneous structures dissolved into the spatially resolved parts and unresolved parts (Agrawal et al., 2001). The resolved parts are modeled using kinetic theory-based TFM, while the unresolved parts are solved using filtered sub-grid constitutive models. The constitutive relations of the filtered drag force and solid stress proposed by Igci et al. (2008) are solved as a function

2

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

of filtered solid volume fraction, Froude number and dimensionless filter length based on the terminal settling velocity of particles (Ozarkar et al., 2015; Sarkar et al., 2016). The filtered drag model proposed by Parmentier et al. (2012) is originated from the filtered relative velocity due to spatial inhomogeneities using the filtered volume averaged TFM and KTGF. The filtered drag force is correlated with the filtered solid volume fraction and dimensionless filter length based on the bed hydraulic diameter. The filtered dragstress model proposed by Scherabueaer et al. (2014) has a very similar form to the method used in the fTFM model (Agrawal et al., 2001), and the drag coefficient and solid viscosity are correlated to dimensionless filter length based on the characteristic length scale. The anisotropic filtered model is expressed as a function of solid volume fraction and dimensionless filter length based on the terminal settling velocity (Cloete et al., 2018). The enhanced filtered model proposed by Gao et al. (2018) is evaluated by applying it to the coarse grid simulations of bubbling fluidized bed, turbulent fluidization, circulating fluidized bed and dilute transport. The filtered stress model has taken the effect of axial gas pressure gradient over the domain on the averaged gas-particles axial flow into account (Mouallen et al., 2019). The gas pressure gradient is enforced through the initial periodic boundaries, indicating domain parameters should be taken into account in the accurate filtered models. The heterogeneous structure-based drag models, including energy minimization multiscale (EMMS) method (Li and Kwauk, 2003), EMMS scheme (Nikolopoulos et al., 2010), EMMS-based model (Shah et al., 2011), EMMS-TFM model (Zeneli et al., 2015), CD-Lab model (Schneiderbauer et al., 2013), structure-dependent drag model (Wenming et al., 2017) and cluster structure dependent (CSD) drag model (Shuai et al., 2014), assume that the gasparticle mixture consists of the heterogeneous structures with the particle-rich dense phase as the form of clusters and the gas-rich dilute phase as the form of dispersed particles. The EMMS method assumes that the suspending and transporting energy is minimum, and is known as a stability condition (Li and Kwauk, 2003). The drag forces are solved by the set of equations with the stability condition. The drag coefficient is expressed by heterogeneity index which is correlated with solid volume fraction only based on Wen and Yu empirical correlation (Wen and Yu, 1966). On the other hand, the empirical correlations of cluster diameter proposed by Gu and Chen (1998) and Harris et al. (2002) are taken as an additional equation. The former empirical correlation is used in EMMS scheme model to determine the heterogeneity index (Nikolopoulos et al., 2010), and the latter empirical equation of cluster diameter is used in the EMMS-based model to determine the correction factor (Shah et al., 2011) and the structure-dependent drag model to calculate the heterogeneous index (Wenming et al., 2017). The analytic equation for cluster diameter proposed by Subbarao (2010) is used in the EMMS-TFM model as a function of cluster volume fraction (Zeneli et al., 2015). Therefore, the simulated drag coefficient of heterogeneous structures depends on empirical correlations or analytic equation of cluster diameters. The CD-Lab model proposed by Schneiderbauer et al. (2013) assumes the exponential distribution of probability density function of cluster diameter. The cluster size is expressed as a function of /ds , where  is the grid spacing, indicating that the determination of heterogeneous structure parameters depends on grid size. The CSD drag model predicts drag coefficient of heterogeneous structures using a set of nonlinear conservation equations (Shuai et al., 2014; Li et al., 2018). The cluster diameter is solved from the equation set without additional empirical or analytical correlations. It is acknowledged that the formation of clusters is mainly caused by the surrounding gas and particles interactions (Goldhirsch and Zanetti, 1993; Wylie and Koch, 20 0 0; Cocco et al., 2010; McMillanet al., 2013) and the inelastic collisions of particles (Luding et al., 1998; Gustavsson et al., 2014). Experi-

mental measurements showed that the heterogeneous structure exhibits temporal correlation with the change of time (Cahyadi et al., 2017). However, the heterogeneous structure- based drag models mentioned above do not take the interactions of collisions of particles and the rates of accumulation of momentum of clusters and dispersed particles with respect to time into account. Hence, the contributions of the interactions of collisions and the rates of momentum accumulation of clusters and dispersed particles on drag forces of heterogeneous structures are necessary. Two kinds of flow structures of gas-particles suspension exist in the computational cell: homogenous flow and heterogeneous flow with clusters, as shown in Fig. 1 in this study. If clusters exist at the state I in the computational cell, the heterogeneous structure-based drag models must be modified to account for the effect of heterogeneous structures on the gas-particles interactions when applying TFM to coarse meshes. On the other hand, when no clusters exist at the state II in the computational cell, and particles are distributed homogenously, the homogenous drag coefficients are predicted by Huilin–Gidaspow drag models (ANSYS, Inc., 2011). Hence, the criterion distinguishing from the heterogeneity of clusters and the homogeneity of particles suspension in the computational cells is essential. On the other hand, the measured particle cluster properties, including cluster diameter, velocity and volume fractions, are reviewed by Harris et al. (2002) and Cahyadi et al. (2017) in CFB risers ranging from laboratory to industrial scale. Therefore, the comparison of experimental data with numerical simulations using the filtered multiscale models and heterogeneous structure-based drag models is needed. In present study, a dynamic cluster structure-dependent (DCSD) drag model is proposed to take the interactions of inter-particle collisions and gas-particles hydrodynamic interactions as well as the temporal accelerations of momentum transports on temporalspatial heterogeneous structure into account. The criterion of the existence of clusters is proposed based on the minimization of energy dissipation rates of the gas-particles interactions and collisional interactions of particles based on KTGF in the computational cells. The intermittency factor is defined to describe intermittent flow of heterogeneous structures. The profiles of cluster solid volume fraction, cluster diameter and intermittency factor are compared to experimental results. The compromise in competition of hydrodynamic interactions and interactions of collisions on clusters existence is analyzed. 2. Euler-Euler gas-solid two-fluid model with KTGF The Euler-Euler TFM is based on the mass and momentum conservation laws of gas phase and solid phase (Gidaspow, 1994). The conservation equations of mass for gas phase (i = g) and solid phase (i = s) are expressed as

∂ ( ε ρ ) + ∇ · ( ε i ρi u i ) = 0 ∂t i i

(1)

The momentum conservation equation of gas phase is expressed as

∂ (ε ρ u ) + ∇ · (εg ρg ug ug ) = −εg ∇ p + εg ∇ · τg + βgs (us − ug ) ∂t g g g + ε g ρg g (2) where τ g is the stress tensors of gas phase. The momentum conservation equation of solid phase is as follows (Gidaspow, 1994)

∂ (ε ρ u ) + ∇ · (εs ρs us us ) = −εs ∇ p − ∇ ps + εs + ∇ · τs ∂t s s s + βgs (ug − us ) + εs ρs g

(3)

where β gs is the drag coefficient. ps is the solid pressure. The solid phase stress tensor τ s is expressed by Eq. (T1.2), as shown

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

3

Fig. 1. Homogenous structure and heterogeneous structure with the dense phase (ε c is less than ε g ) and the dilute phase (ε d is larger than ε g ) in the computational cell.

in Table 1, where the particles shear viscosity μs and bulk viscosity ξ s are modeled as a function of granular temperature θ using KTGF. The conservation of particles fluctuating energy for granular temperature is expressed as (Gidaspow, 1994)

3 2





∂ (ε ρ θ ) + ∇ · (εs ρs θ )us = (−∇ ps I + τs ) : ∇ us ∂t s s +∇ · (ks ∇θ ) − s + φs + Dgs

Table 1 Constitutive equations of gas-particles two-fluid model. (a) Gas phase stress

   1 τg = μ g ∇ u g + ( ∇ u g ) T − ( ∇ · u g ) I 3

(T1.1)

(b) Solid phase stress

(4)

Hence, the computational cell parameters of gas volume fraction ε g , gas pressure p, granular temperature θ , and velocities of gas phase ug and solid phase us are solved from TFM with KTGF in the computational cells. 3. Dynamic cluster structure-dependent drag model The formation of heterogeneous structure depends on gasparticles interactions which relate to slip velocity between gas velocity and solid velocity. For fluid catalytic cracking (FCC) particles, the solid velocity is close to the gas velocity in the direction of the flow (Grace et al., 1997; Gidaspow, 1994; Huilin, 2017). Thus, onedimensional flow of heterogeneous structures is assumed along the mean flow direction (vertical z direction) in a first approximation. Although one-dimensional equations of heterogeneous structures are easier to obtain from the momentum equations of gas and solid phases, the equations are solved to be useful for further investigations of multidimensional equations of hydrodynamics of heterogeneous structures. The heterogeneous structure at the state I, as shown in Fig. 1, is dissolved into the dilute phase where the gas volume fraction is larger than ε g , dense phase in which the gas volume fraction is less than ε g and the interface between the dilute phase and the dense phase in the computational cell volume V. For the state I, the dense phase volume part Vc consists of clusters with a mean cluster diameter dc , gas superficial velocity Ug,c , clusters superficial velocity Us,c , gas volume Vg,c and granular temperature of particles within the clusters θ c . The dilute phase volume part Vd of dispersed particles is described by the gas superficial velocity Ug,d , gas volume Vg,d , dispersed particles superficial velocity Us,d and dispersed particles granular temperature θ d . The dense phase volume fraction of the computational cell defines f=Vc /V. The gas volume fractions of the dense phase and the dilute phase are ε c =Vg,c /Vc and ε d =Vg,d /Vd in the computational cell.

  1 τs = μ s { ∇ u s + ( ∇ u s ) T − ( ∇ · u s ) I } + ξ s ∇ · u s 3

(T1.2)

(c) Solid pressure

ps = εs ρs θ [1 + 2g0 εs (1 + e )]

(T1.3)

(d) Shear viscosity of solids



√

θ 10ρs ds πθ + π 96(1 + e )εs g0

2 4 × 1 + go εs ( 1 + e ) 4 2 ε ρs ds g0 (1 + e) 5 s

μs =

(T1.4)

5

(e) Bulk solids viscosity

ξs =

4 2 ε ρs ds g0 (1 + e ) 3 s



θ π

(T1.5)

(f) Thermal conductivity of particles

ks =

√

2 25ρs ds πθ 6 1 + ( 1 + e )g0 εs 64(1 + e )g0 5 +2εs2 ρs ds g0 (1 + e )

 1/2 θ π

(T1.6)

(g) Dissipation fluctuating energy

 γs = 3 1 − e2 εs2 ρs g0 θ

 4 ds



θ − ∇ · us π

 (T1.7)

(h) Radial distribution function

 g0 = 1 −

 ε 1/3 −1 s εs max

(T1.8)

(i) Drag coefficient

βgs = ϕgs =

 βDCSD βHG

BEV Ndf = minimum BEV Ndf = minimum

arctan[150 × 1.75(0.2 − εs )]

π

+ 0.5

(T1.9)

(T1.10)

(j) Rate of energy dissipation per unit volume

3.1. Momentum equations of dilute phase and dense phase Rather than solving the complete transport equations of dispersed particles in the dilute phase, we assume that the offdiagonal elements of the gradient of the dispersed particles stress tensor are zero (Gidaspow, 1994). The one-dimensional transient momentum equation of dispersed particles is expressed by

Dgs =

ds ρs √

4 πθ go



18μg ds2 ρs

2 |ug − us |2

(T1.11)

(k) Exchange of fluctuating energy between gas and particles

φs = −3βgs θ

(T1.12)

4

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

Eq. (T2.1), as shown in Table 2 (Huilin, 2017; Li et al., 2018). The second term on the right hand side is the drag force in the dilute phase. The forth term is the dispersed particles pressure, and the last term is momentum exchange between the dilute phase and the dense phase. Because the density of particles is greater than that of gas phase, Eq. (T2.1) can be rewritten to

nd Fd = (1 − f )(1 − εd )(ρs − ρg ) g + asd,l + asd,c



 Us,d ∂ p ∂ ps,d Us,c + ( 1 − f ) ( 1 − εd ) + + βdc − ∂z ∂z 1 − εc 1 − εd

where asd,l and asd,c are the temporal acceleration and convective acceleration of dispersed particles, and they can be deduced from the given equations. The temporal acceleration and convective acceleration represent the rate of velocity change with respect to time and the spatial gradients of velocity of dispersed particles in the dilute phase. β dc is the drag coefficient of the solids-solids momentum transfer between the dilute phase and the dense phase. An equation proposed by Syamlal et al. (1993) is used to model drag coefficient between the dilute and dense phases. We assume that the off-diagonal elements of the gradient of the solid stress tensor of clusters are zero in a first approximation. The one-dimensional transient momentum equation of clusters is expressed by Eq. (T2.2), as shown in Table 2 (Huilin, 2017; Li et al., 2018). The second and third terms on the right hand side represent the drag forces in the dense phase and the interface between the dilute phase and the dense phase. Because the gas density is less than that of particles, Eq. (T2.2) can be rewritten to



 Us,d Us,c ∂ p ∂ ps,c + f ( 1 − εc ) + − βdc − ∂z ∂z 1 − εd 1 − εc

(6)

where asc,l and asc,c are the cluster temporal acceleration and convective acceleration in the dense phase, and represent the rate of cluster velocity change with respect to time and the spatial gradients of velocity of clusters due to the change of position. Both accelerations can be deduced from the given equations. Similar to the assumptions used in the momentum equations of dispersed particles and clusters, the one-dimensional transient momentum equations of gas phase in the dilute phase and the dense phase along the direction of the flow are expressed by Eq. (T2.3) and (T2.4) (Huilin, 2017; Li et al., 2018). Combining with them, the gas momentum equation becomes as follows:

  nc Fc nd Fd ni Fi = + + ρg (agd,c − agc,c ) + (agd,l − agc,l ) f εc ( 1 − f )εd ( 1 − f )εd (7) where the gas temporal accelerations agc,l and agd,l , and gas convective accelerations agc,c and agd,c of the dense phase and dilute phase can be deduced from the given equations. 3.2. Mass balance of gas phase and solid phase The gas and particles mass balance equations in the dense phase and dilute phase are expressed in the computational cell as following.

f ρgUg,c + (1 − f )ρgUg,d = εg ρg ug

(8)

f ρsUs,c + (1 − f )ρsUs,d = (1 − εg )ρs us

(9)

The gas volume fraction requires that the sum of volume fraction of gas and particles be unity in the computational cell, that is,

εg = f εc + ( 1 − f )εd

The algebraic expression for granular temperature proposed by Syamlal et al. (1993) assumes that the granular energy is dissipated locally; neglecting the convection and diffusion contributions, and retaining the generation and dissipation terms. In present study, the algebraic expression model is extended to predict granular temperatures of dispersed particles θ d and granular temperature of particles within the clusters θ c in a first approximation.

 (5)

nc Fc + ni Fi = f (1 − εc )(ρs − ρg ) g + asc,l + asc,c

3.3. Fluctuating energy equations for dispersed particles and clusters

(10)

−K1,d ∂ 2K4,d ∂ z

θd =



+



K12,d

(1 − f )Us,d ( 1 − εd )

K2,d + K5,d

4K42,d





⎫    2 ⎬2 2K3,d ∂ (1 − f )Us,d ∂ 2 (1 − f )Us,d + K6,d ∂ z ( 1 − εd ) ∂ z 2 ( 1 − εd ) ⎭ (11)

 θc =

−K1,c ∂ 2K4,c ∂ z



+



K12,c 4K42,c



f Us,c

( 1 − εc ) K2,c + K5,c



⎫    2 ⎬2 f Us,c 2K3,c ∂ f Us,c ∂2 + K6,c ∂ z (1 − εc ) ⎭ ∂ z 2 ( 1 − εc ) (12)

where the coefficients Ki (i = 1–6) are given in Table 2. e is the coefficient of restitution of particles. 3.4. Equations of energy balance and minimization of energy dissipation rate The kinetic energy equation of solid phase was derived according to the kinetic-theory-based model (Gidaspow, 1994; Fox, 2014). For flow of heterogeneous structures, the granular energy Es consists of the energy component of dispersed particles in the dilute phase Es,d and the energy component of clusters in the dense phase Es,c according to the velocities of the dilute phase us,d and dense phase us,c .

Es,d =

1 3 u ·u + θ 2 s,d s,d 2 d

(13a)

Es,c =

1 3 us,c · us,c + θc 2 2

(13b)

The conservation energy equation of solid phase is expressed by Eq. (T2.5), as shown in Table 2. The fifth term on the left hand side represents the granular temperature flux components of the dilute and dense phases. A similar operation is performed to analyze the gas energy Eg in the dilute phase and dense phase. The conservation energy equation of gas phase is expressed by Eq. (T2.6). The fourth term on the left hand side represents gas heat fluxes of the dilute and dense phase. Then, the addition of these equations yields the energy balance of gas-particles suspension as follows:

 ∂ ∂  2 2 f (ρsUs,c Es,c + ρgUg,c Eg,c ) + (1 − f ) ρsUs,d Es,d + ρgUg,d Eg,d ( E + Eg ) + ∂t s ∂z     (1 − f )Us,d ∂p ∂ f Us,c ∂ + ( εg ug + εs us ) + ps,c + ps,d + q + qs,c ∂z ∂z 1 − εc 1 − εd ∂ z s,d      ∂ 2 + q + qg,c + f 2 ρsUs,c + (1 − f ) ρsUs,d g + f ρgUg,c + (1 − f )ρgUg,d g ∂ z g,d

 = [nd Fd Ud + nc FcUc + ni (1 − f )FiUi ] + ( s,c + s,c ) + g,d + g,c (14)

Es = f (1 − εc )ρs Es,c + (1 − f )(1 − εd )ρs Es,d

(15a)

Eg = f εc ρg Eg,c + (1 − f )εd ρg Eg,d

(15b)

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

5

Table 2 Model equations of DCSD drag coefficient. (1) One-dimensional transient momentum equation of dispersed particles

   ∂ (1 − f )Us,d ∂ (1 − f )2 ρsUs,d + (1 − f )2 ρsUs,d ∂t ∂z ( 1 − εd ) = −(1 − f )(1 − εd )

 ∂ ps,d Us,d ∂p Us,c + nd Fd − (1 − f )(1 − εd )ρs g − − βdc − ∂z ∂z 1 − εc 1 − εd

(T2.1)

(2) One-dimensional transient momentum equation of clusters

  ∂ 2 ∂ f Us,c ( f ρsUs,c ) + ( f 2 ρsUs,c ) ∂t ∂z ( 1 − εc ) = − f ( 1 − εc )

 Us,d ∂p ∂ ps,c Us,c + nc Fc + ni Fi − − f (1 − εc )ρs g + βdc − ∂z ∂z 1 − εd 1 − εc

(T2.2)

(3) One-dimensional transient momentum equations of gas in the dilute and dense phases

   ∂ (1 − f )Ug,d ∂ (1 − f )2 ρgUg,d + (1 − f )2 ρgUg,d ∂t ∂z εd ∂p = − ( 1 − f )εd − nd Fd − ni Fi − (1 − f )εd ρg g ∂z   fU ∂ 2 ∂ ( f ρgUg,c ) + ( f 2 ρgUg,c ) g,c ∂t ∂z εc ∂p = − f εc − nc Fc − f εc ρg g ∂z

(T2.3)

(T2.4)

(4) One-dimensional kinetic energy equations of particles and gas

 ∂ Es ∂  2 ∂p + f ρs Es,cUs,c + (1 − f )2 ρs Es,d Us,d + εs us ∂t ∂ z ∂z      (1 − f )Us,d ∂ f Us,c ∂ + ps,c + ps,d + (q + qs,c ) + f 2 ρsUs,c + (1 − f )2 ρsUs,d g ∂z 1 − εc 1 − εd ∂ z s,d  

 (1 − f )Us,d f Us,c = nd Fd + nc Fc + ni (1 − f )FiUi + s,d + s,c 1 − εd 1 − εc

(T2.5)

 ∂ Eg ∂  2 2 + f ρg Eg,cUg,c + (1 − f ) ρg Eg,d Ug,d ∂t ∂z   ∂p ∂ + εg ug + (qg,d + qg,c ) + f ρgUg,c + (1 − f )ρgUg,d g ∂ z ∂ z   (1 − f )Ug,d f Ug,c = − nd Fd + nc Fc + ( g,d + g,c ) εd εc

(T2.6)

(5) Drag forces of gas and particles in the dilute and dense phases

    Ug,d Us,d  Ug,d Us,d π ds2 ρg 24 3.6  + − − 1 − εd  εd 1 − εd 8εd2.65 Red Re0d.313  εd

Fd =

Fc =

Fi =

(T2.7)

    Ug,c π ds2 ρg 24 3.6 Us,c  Ug,c Us,c + 0.313  − − 2.65  Rec εc 1 − εc εc 1 − εc 8εc Rec

(T2.8)

     Ug,d  Ug,d π dc2 ρg εd2 24 3.6 Us,c Us,c   + − ( 1 − f ) − 1 − f ( ) 2.65  Rei 1 − εc εd 1 − εc Re0i .313  εd 8 (1 − f )

(T2.9)

(d) Drag coefficient of particles between the dilute phase and the dense phase



βdc =





Us,d  3(1 + e )ρs (1 − εd )(1 − εc ) 1 3(2 − εd − εc )  Us,c +  1 − εc − 1 − εd  ds εg 2εg2

(T2.10)

(6) Solid pressure of the dilute phase and dense phase





ps,d = ρs (1 − f )(1 − εd ) 1 + 2g0,d (1 − f )(1 − εd )(1 + e ) θd

(T2.11)

ps,c = ρs f (1 − εc )[1 + 2g0,c f (1 − εc )(1 + e )]θc

(T2.12)

(7) Number density of cluster, dispersed particles and interface

nc =

f ( 1 − εc ) , π ds3 /6

nd =

(1 − f )(1 − εd ) , π ds3 /6

ni =

f

π dc3 /6

(T2.13)

(8) Radial distribution function of dilute and dense phases

g0,d =

3 3 εs,1/max εs,1/max , g0,c = 1/3 1/3 3 {εs,1/max − [(1 − f )(1 − εd )] } {εs,max − [ f (1 − εc )]1/3 }

(T2.14)

(continued on next page)

6

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126 Table 2 (continued) (9) Collisional energy dissipations of dispersed particles and clusters 1/3 3(1 − e2 )[(1 − f )(1 − εd )] ρs εs,max θd 2

s,d =

4 ds

1/3 1/3 εs,max − [(1 − f )(1 − εd )]

1/3 3(1 − e2 )[ f (1 − εc )] ρs εs,max θc 2

s,c =





1/3 1/3 {εs,max − [ f ( 1 − εc )] }

4 ds





θd ∂ (1 − f )Us,d − π ∂ z 1 − εd

θc ∂ f Us,c − π ∂ z 1 − εc

 (T2.15)

 (T2.16)

(10) Rate of energy dissipations of gas in the dilute and dense phases



g,d =

4 ∂ ( 1 − f ) ε d μg 3 ∂z

g,c =

4 ∂ f ε c μg 3 ∂z







f Ug,c

(1 − f )Ug,d εd

2 (T2.17)

2 (T2.18)

εc

(11) Coefficients

K1,d = 2(1 + e )ρs g0,d ,

K2,c =

K2,d =

K3,c =

K3,d

K1,c = 2(1 + e )ρs g0,c

8ds ρs (1 + e ) f (1 − εc )g0,c

√  3 π   √ π 8 f (1 − εc )g0,c (1 + e ) ds ρs − √ [1 + 0.4(1 + e )(3e − 1 ) f (1 − εc )g0,c ] + 6 3 (3 − e ) 5 π

(T2.20)

8ds ρs (1 + e )(1 − f )(1 − εd )g0,d

√    √3 π  8(1 − f )(1 − εd )g0,d (1 + e) π  ds ρs − 1 + 0.4(1 + e )(3e − 1 )(1 − f )(1 − εd )g0,d + √ 6 3 (3 − e ) 5 π

(T2.21)

√

ds ρs 4

ds ρs = 4

π [1 + 0.4(1 + e)(3e − 1) f (1 − εc )g0,c ] 8 f (1 − εc )g0,c (1 + e) + √ 3 (3 − e ) 5 π



12 1 − e2 ρs g0,c , √ ds π

K5,c =

12 1 − e2 f (1 − εc )ρs g0,c , √ ds π



 (T2.22)

√    π 1 + 0.4(1 + e)(3e − 1)(1 − f )(1 − εd )g0,d 8(1 − f )(1 − εd )g0,d (1 + e ) + √ 3 (3 − e ) 5 π

K4,c =

K6,c =

(T2.19)

K4,d =



12 1 − e2 ρs g0,d √ ds π





96 1 − e2 f (1 − εc )ρs g0,c √ ds π

K5,d =

(T2.24)



12 1 − e2 (1 − f )(1 − εd )ρs g0,d √ ds π



2 ,

(T2.23)

K6,d =



96 1 − e2 (1 − f )(1 − εd )ρs g0,d √ ds π

In this equation, the gas energy components of the dilute phase Eg,d and the dense phase Eg,c are unclosed, and must be found from the transport equations, leading to Eq. (14) is more complicated than in the gas energy conservation equation. But, we observe that the terms at the right hand side represent the rate of energy dissipations of gas-particles suspension through three contributions by hydrodynamic forces, interactions of collisions of particles and gas viscous forces of heterogeneous structures.

   = [nd FdUd + nc FcUc + ni FiUi (1 − f )] + s,d + s,d + g,d + g,c (16) To keep the energy dissipation rate as low as possible, particles tend to aggregate into clusters to reduce the gas flow resistance and the collisional energy dissipations, while gas tends to bypass the clusters instead of penetrating through clusters (Gidaspow, 1994; Huilin, 2017). Hence, the criterion of the existence of heterogeneous structures is the energy dissipation rate per unit mass

(T2.25)

2 (T2.26)

of mixture is minimum:

Ndf =

 = dr + co + vi = minimum ε g ρg + ( 1 − ε g ) ρs

(17)

This is an extreme value equation of the rate of energy dissipation including the energy dissipation component of drag forces dr , the collisional energy dissipation component co and gas viscous energy dissipation component vi . It is important to keep in mind that we solve an extreme value Eq. (17) of the minimization of energy dissipation rate, instead of solving the balance Eqs. (14) of granular energy Es and gas energy Eg to determine the heterogeneous structure parameters of the dense phase and dilute phase in the computational cell. Like energy balance Eq. (14), the extreme value Eq. (17) takes as a closure equation for heterogeneous structures.

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

3.5. Conditioned extreme value equation based on bivariate extreme value (BEV) theory For flow of heterogeneous structures with dispersed particles in the dilute phase and clusters in the dense phase, there are ten independent variables at the state I, as shown in Fig. 1. On the other hand, we have a set of equations with an extreme value Eq. (17) and eight balance equations including three momentum Eqs. (5)–(7), three continuity Eqs. (8)–(10) and two granular temperature Eqs. (11) and (12). Hence, the set of transient nonlinear equations is unclosed to solve ten heterogeneous structure parameters in the computational cell. The dependence function with two independent variables is described by the bivariate extreme value (BEV) theory (Kotz and Nadarajah, 20 0 0). The BEV method determines an extreme point of the dependence function with the change of two independent variables if it has a solution. Thus, the values of two independent variables are determined from the minimum value point of the dependence function. In present study, among ten heterogeneous structure parameters the gas volume fractions of the dense phase εc and the dilute phase ε d are chosen as two independent variables of the dependence function Ndf . Using BEV theory, the dependence function Ndf from Eq. (17) is expressed as follows ε , p, u , u , θ

BEV Ndf = (dr + co + vi )|εg ∈R1 , εg c ∈Rs 2 = minimum d

(18)

This is a conditioned extreme value equation, and it gives the relation between the minimum value point of the dependence function Ndf and two heterogeneous structure parameters ε d and εc in the ranges of the dilute phase R1 and the dense phase R2 according to five computational cell parameters in the computational cell. The ranges of R1 and R2 are R1 ∈[ε g , 1.0] and R2 ∈[ε mf , ε g ], where ε mf is the gas volume fraction at the minimum fluidization. The conditioned extreme value Eq. (18) indicates that the dependence function Ndf has a minimum value point in the range of R1 and R2 if the heterogeneous structure exists in the computational cell. Thus, the heterogeneous structure parameters ε d and ε c of the computational cell are determined at the minimum value point. Furthermore, once the existence of clusters is identified in the computational cell, the remaining eight heterogeneous structure parameters (i.e., dc , f, Ug,d , Us,d , θ d , Ug,c , Us,c and θ c ) are obtained by solving eight balance equations. Hence, the set of transient nonlinear equations is determined. On the other hand, Eq. (17) has no a minimum value point, and the rate of energy dissipation Ω is monotonically increased with the change of gas volume fractions of the dilute phase ε d in the range R1 and the dense phase ε c in the range R2 . This identifies a homogenous flow (HF) of gas and particles mixture at the state II, as shown in Fig. 1 in the computational cell. ε , p, u , u , θ

BEV Ndf = (dr + co + vi )|εg ∈R1 , εg c ∈Rs 2 = minimum d

(19)

The solutions of Eqs. (18) and (19) determine the heterogeneity of clusters and the homogeneity of particles suspension in computational cells. They indicate that the heterogeneous structures exist when the dependence function of energy dissipation rate has a minimum value point. The existence of the minimum value point of the energy dissipation rate depends on the computational cell parameters, suggesting that the direction of evolution of heterogeneous structure minimizes its energy dissipation rate compatible with the constrains of computational cell parameters. 3.6. Dynamic cluster structure-dependent (DCSD) drag model From gas momentum balance equation with no gas acceleration and gravity in the dense phase and dilute phase (Huilin, 2017), the

7

DCSD drag coefficient is expressed as follows:

βDCSD =

nc Fc + nd Fd + ni Fi |ug − us |

BEV Ndf = minimum

(20)

On the other hand, the homogeneous drag coefficient is calculated using Huilin–Gidaspow drag correlation with the switch function ϕ gs :



βHG = ϕgs

150εs (1 − εg )μg

+(1 − ϕgs )

1.75εs ρg |us − ug | + ds

εg ds 2 3CD εs εg ρg |us − ug | 4ds



εg −2.65 NdfBEV = minimum (21)

From Eqs. (20) and (21), the drag coefficient in the computational cell is modeled as follows:



βgs =

βDCSD βHG

BEV Ndf = minimum BEV Ndf = minimum

(22)

3.7. Implementation of DCSD drag model in MFIX In this study of numerical simulations, the algorithm is based on the modified CFD MFIX- TFM which is fully described by NETL (National Energy Technology Laboratory, DOE, USA) (Syamlal et al., 1993). In the MFIX-TFM code, the set of equations of the DCSD drag coefficient are implemented in the drag_gs subroutine, located in drag_gs.f. The discretization of the momentum Eqs. (5)–(7) is performed using the first-order scheme. The discretized formulas are solved using equation solvers which are available in MFIXTFM code. Two solvers are used. One is the linear equation solver. Another is the equation solver for the solution of minimization energy dissipation rate using the steepest descent method. We use the iterative solvers as we have an initial guess for the solution. An optimum degree of convergence has been determined and is controlled by the specified iterations inside the equation solvers given in the mfix.dat file. To predict flow behavior of gas and particles in CFB riser in this study, the experimental results measured by Wei et al. (1998) are tested. The riser diameter and height are 0.186 m and 8.0 m. The particles diameter and density are 54 μm and 1398 kg/m3 . The coefficient of restitution of particles is 0.95 in numerical simulations. The gas is supplied at the bottom of the riser. The gas and solid pass through the outlets at the top of the riser. The exit is set to continuous outflow boundary conditions for both gas and solid phases. The exit gas pressure is assumed to be 101,330 N/m2 . At the wall, a no-slip condition is applied for the gas phase. The boundary conditions of Johnson and Jackson (1987) for the tangential velocity and granular temperature are applied. Initially, the velocities of gas and solid phases are zero in the riser. 4. Simulations and discussions 4.1. Comparison of simulations and experiments The numerical simulations by means of the increments of grid numbers in both the radial and axial directions are performed using DCSD drag model and Gidaspow drag model (Gidaspow, 1994). Fig. 2 shows the profiles of computed time-averaged solid volume fractions along radial direction. The computational domains consist of 16 × 216 and 20 × 284 grids (radial direction × axial direction). Both the grid number cases predict similar distributions with the low solid volume fractions at the center and high solid volume fractions near the wall. The solid volume fraction error percentages of the coarse and fine grid simulations as compared

8

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

Fig. 2. Distributions of solids volume fractions along radial direction.

with experiments are 4.06% and 3.34% using DCSD drag model, respectively. This indicates that both the grid number cases are sufficiently fine for providing reasonably grid independence results. On the other hand, the solid volume fraction error percentage of the finer grids case as compared with experiments is 30.31% using Gidaspow drag model. With the decrease grid number, the error percentage increases since Gidaspow drag model does not account for any heterogeneity of gas-particles suspension. The continuous increase in mesh density may lead to better results in comparison with experiments. However, the computational time is still a significant restriction when using a finer grid. In this study, the computational domain consisting of 16 non-uniform grids in radial direction and 216 uniform grids in axial direction is selected to use. The simulated instantaneous solids volume fractions are shown in Fig. 3 using DCSD drag model and Gidaspow drag model (Gidaspow, 1994) at the inlet gas superficial velocity and solid mass flux of 3.25 m/s and 98.8 kg/m2 s. Both cases show that the characteristic feature of the flow is the oscillating motion

Fig. 4. Profile of gas and solid volume fractions along riser height.

of particles with a high solids volume fraction near the walls and low at the center of riser. It further illustrates that the solids volume fraction is high near the bottom due to the continual collection of particles. These temporal variations of solid volume fraction are perceived clearly as the formation of clusters near the walls comparing to the clusters existence at the center. The time averaged solid volume fractions are calculated from 20 to 60 s, and the comparison of simulated time-averaged axial profiles of gas volume fraction and experimental data is shown in Fig. 4. All simulations and experiments are low at the bottom and high at the top region. It further illustrates that the simulated and measured solid volume fractions are lower at the center than that near the wall. Using DCSD drag model and experimental data, there is a transition from low to high and high to low at a height of 0.5–2 m from the inlet, and the difference between the simulations and experiments may be mostly attributed to that a real 3-D inlet geometry described by a 2-D uniform inlet condition.

Fig. 3. Instantaneous solids volume fractions using (a) DCSD drag model and (b) Gidaspow drag model at the superficial gas velocity of 3.25 m/s.

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

9

Fig. 5. Profile of minimization energy dissipations of drags of the co-current up-flow of gas and particles at the riser height and time of z = 4.5 m and t = 35 s.

4.2. Distribution of minimization of energy dissipation rate of heterogeneous structure The criterion of cluster existence indicates that the energy dissipation rate is minimum as the gas volume fractions of ε d and εc vary in the ranges of R1 for the dilute phase and R2 for the dense phase. The profile of instantaneous energy dissipation rate is shown in Fig. 5 at the computational cell (a) at ug =2.87 m/s and us =1.23 m/s near the center (x/D = 0.089) and (b) at ug =0.295 m/s and us =0.175 m/s near the wall (x/D = 0.45). Both the computational cell parameters of instantaneous gas velocity and solid velocity are positive, indicating a co-current up-flow of gas and particles. In the computational cell (a), several extreme points of energy dissipation rates are found with the change of ε c and ε d . The multi-extreme point behavior originates from the contributions of hydrodynamic interactions and collisional collisions of particles. The minimum value point is found from these extreme points, indicating the clusters exist in the computational cell. The values of ε d and ε c at the minimum value point are 0.991 and 0.860. The remaining eight heterogeneous structure parameters are calculated by means of the equation set, and are listed in Table 3. On the other hand, for the computational cell (b) the energy dissipation rate is monotonically increased with the increase of ε c and decrease of ε d . There are no extreme points from Eq. (18) with the change of ε d in the range R1 and ε c in the range R2 , indicating that the HF of gas-particle suspension exists in the computational cell.

Fig. 6 shows the distributions of instantaneous energy dissipation rate at the computational cell (a) at ug =0.41 m/s, us = −0.32 m/s and x/D = 0.089 and (b) at ug =1.103 m/s, us =−1.329 m/s and x/D = 0.12. Both cases have a positive gas velocity and negative solid velocity, indicating a counter-current flow with gas upflow and particles down-flow in the computational cells. In the computational cell (a), the distribution of energy dissipation rate has several extreme points with the change of ε c and ε d . One minimum value point is determined, and the values of gas volume fractions of dilute phase and dense phase are 0.982 and 0.801 at the minimum value point. The gas superficial velocities of the dilute phase and dense phase shown in Table 3 are positive. The dispersed particles superficial velocity of the dilute phase Us,d is positive and the cluster superficial velocity of the dense phase Us,c is negative, suggesting the clusters flow downward and the dispersed particles flow upward in the computational cell. On the other hand, the energy dissipation rate in the computational cell (b) increases with the increase of ε c and ε d , indicating the HF of gas and particles suspension in the computational cell. The profiles of instantaneous energy dissipation rate are shown in Fig. 7 at the computational cell (a) at ug =−0.27 m/s, us = −0.5 m/s and x/D = 0.19 and (b) at ug =−0.142 m/s, us =−0.387 m/s and x/D = 0.38. The computational cell parameters of gas velocity and solid velocity are negative, indicating a co-current down-flow of gas and particles. In the computational cell (a), the heterogeneous flow is identified from the extreme points of energy dissipa-

Table 3 Computational cell parameters from TFM and heterogeneous structure parameters at the minimum point from DCSD drag model. Parameters Computational cell parameters

Heterogeneous structure parameters

At minimum point Equation set

ug (m/s) us (m/s)

εg θ (m/s) 2 p/z(Pa/m) εd εc f dc (mm) Us,d (m/s) Us,c (m/s) Ug,d (m/s) Ug,c (m/s) θ d (m/s)2 θ c (m/s)2

Co-current up-flow

Counter-current flow

Co-current down-flow

2.87 1.23 0.892 0.00735 −230.9 0.991 0.860 0.755 4.32 1.11 × 10−3 0.175 3.78 2.17 1.61 × 10−1 4.24 × 10−4

0.41 −0.32 0.817 0.0031 −134.6 0.982 0.801 0.912 1.76 7.33 × 10−2 −7.09 × 10−2 2.69 0.11 1.87 × 10−3 5.26 × 10−4

−0.34 −0.43 0.785 0.00051 −361.0 0.964 0.385 0.309 0.684 −5.27 × 10−3 −0.288 −0.256 −0.291 6.91 × 10−2 4.64 × 10−4

10

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

Fig. 6. Profile of energy dissipations of drags of the counter-current flow of gas and particles at the riser height and time of z = 4.2 m and t = 38 s.

Fig. 7. Profile of minimization energy dissipations of drags of the co-current down-flow of gas and particles at the riser height and time of z = 4.13 m and t = 33 s.

Fig. 8. Instantaneous cluster diameter and solid volume fraction at two computational cells.

tion rate. The gas volume fractions of the dense phase and dilute phase are 0.385 and 0.964 at the minimum value point of energy dissipation rate. The superficial velocities of gas, dispersed particles and clusters in the dense phase and dilute phase, as shown in Table 3, are negative. On the other hand, the HF in the computational cell (b) is identified because no extreme points exist with the change of ε d in the range R1 and ε c in the range R2 in the computational cell

4.3. Predictions of heterogeneous structure parameters The instantaneous cluster diameters dc and solid volume fractions ε s are shown in Fig. 8 at two computational cells. More clusters are predicted near the wall comparing to the formation of clusters in the center of the riser. The riser wall stimulates cluster formation with high solid volume fractions. Simulations show the cluster diameter is larger near the wall than that in the center,

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

11

Fig. 9. Instantaneous gas volume fraction and solids volume fraction of dense phase at two computational cells.

indicating the high solid volume fraction produces denser clusters. When the solid volume fraction is low, the clusters are vanished, indicating the cluster appearance intervals with time. The occurrence numbers of clusters Nc are 53 and 156 at two computational cells. Soong et al. (1993) proposed the guidelines to define clusters, and the value of the local instantaneous solid volume fraction for a cluster is higher than its time-averaged value ε s by a constant δ times the standard deviation (SD), i.e. ε s +δ × SD, where the constant δ is 2.0. The cluster numbers Ndc are determined, and they are 28 and 49, which are less than Nc . On the other hand, the value of constant n was 1.0–1.4 (Manyele et al., 2002). When the constant δ reduces to 1.0, more clusters are detected from simulated instantaneous solid volume fractions in the computational cell. The cluster number Ndc using the constant δ of 1.0 are 67 and 193 which are close to the occurrence number of clusters Nc at two computational cells. The instantaneous clusters solid volume fraction ε s,c of the dense phase and gas volume fraction ε g are shown in Fig. 9 at two computational cells. Simulations show that both clusters solid volume fractions and gas volume fractions vary with time and position. High gas volume fractions lead to fewer clusters. The mean solid volume fractions and standard deviations of clusters are 0.265 and 0.051 near the wall and 0.175 and 0.104 at the center, respectively. These show that the clusters solid volume fraction is larger near the wall than that at the center. Contrary, the standard deviation is large at the center, and low near the wall, indicating the wall suppresses the motion of clusters. The instantaneous cluster velocities, shown in Fig. 10, indicate that the mean cluster velocity and standard deviation

Fig. 10. Simulated instantaneous cluster velocity at two computational cells.

are 0.944 m/s and 0.708 m/s at the center and −0.235 m/s and 0.112 m/s near the riser wall. The simulated negative and positive cluster velocities suggest that clusters travel downward in the vicinity of the riser wall, and flow upward in the center. The instantaneous convective accelerations of clusters and dispersed particles are shown in Fig. 11 at two computational cells. The positive and negative convective accelerations indicate that the clusters and dispersed particles are accelerated or decelerated with time. The accelerations of clusters are significantly higher than that

Fig. 11. Instantaneous convective accelerations of clusters and dispersed particles.

12

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

Fig. 13. Profile of simulated intermittency factors and measured cluster residence time fractions.

Fig. 12. Distribution of intermittency factor at four riser heights.

of dispersed particles because of the large mass of clusters. The variation of convective acceleration causes the significant velocity evolutions through acceleration and/or deceleration of dispersed particles and clusters. 4.4. Intermittency factor and intermittent behaviors of clusters For a fluid flow of laminar-to-turbulent transition, the intermittency defines the time fraction of the flow is turbulent. The clusters occur at the time τ i , and dissolve at the time τ i + 1 . The cluster existence time is τ i , as shown in Fig. 9. The intermittency factor defines the ratio of the total cluster duration time to the total sampling time τ in the computational cell.

ξ=



i=1

τ

τ i

(23)

The intermittency factor measures the existence time fraction of clusters, and reveals the alternation of heterogeneous flow and homogeneous flow in the computational cell. From the definition of intermittency factor, it can be seen that the intermittency factor of 0.0 represents HF of gas and particles suspension and no clusters occur with time. On the other hand, the value of 1.0 indicates the fully heterogeneous structures (FHS) with clusters. The distribution of intermittency factor is shown in Fig. 12 at four heights. The intermittency factor is low at the center, and large toward the wall. The values of intermittency factor are in the range of 0.068–0.105 at the center, and 0.186–0.418 near the wall. Furthermore, the intermittency factor is large at the bottom, and small at the upper. The mean intermittency factor is from 0.094 to 0.246, indicating the flow of heterogeneous structure is intermittent in the riser. From experimental measurements, the residence time fraction of clusters was defined (Wu et al., 1991; Wei et al., 1995; Sharma et al., 20 0 0), and it represents the ratio of the total duration time of clusters to the total sampling time. From the definition, we see that the value of the residence time fraction is equal to intermittency factor using Eq. (23). The simulated intermittency factors and measured cluster residence time fractions are shown in Fig. 13 with the change of solid volume fractions. Both simulated intermittency factors and measured cluster residence time fractions increase with the increase of solid volume fractions. The intermittency index was defined by Brereton and Grace (1993) from measured solid volume fraction fluctuating signals. The intermittency index defines the ratio of the standard deviation of solid volume fraction fluctuation to the standard deviation of solid volume fraction fluctuation for an ideal flow of clusters. Fig. 13 illustrates that the experimental measurements

Fig. 14. Distribution of frequency of cluster occurrence at four riser heights.

from Lints et al. (1993) and Wu et al. (1991) are much high in comparison with experimental results measured by Wei et al. (1995) and Sharma et al. (20 0 0). The measurements from Sharma et al. (20 0 0) appear to be relatively constant and insensitive to changes of solid volume fractions. The discrepancy of numerical simulations and experimental data is most likely due to the criterion used in the identification of clusters. The cluster frequency defines the number of clusters observed per second over the total observation time in the computational cell. The calculated cluster frequency is shown in Fig. 14 at four riser heights from simulated instantaneous cluster solid volume fraction. The cluster frequency is low at the center, and large toward the wall. The values of cluster frequency are in the range of 6.93–9.33 Hz at the center, and 17.14–42.13 Hz near the wall. The measured cluster frequency by Sharma et al. (20 0 0) was about 6.0– 10.0 Hz. These numerical simulations and experiments suggest the occurrence of clusters is a transient phenomenon in the riser. 4.5. Compromise of energy dissipations of gas-particles interactions and inter-particle interactions The formation of heterogeneous structure contributes to the interactions of particles with other particles and gas phase. The inelastic collisions of particles dissipate kinetic energy co , and the hydrodynamic interactions between the gas phase and solid phase

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

13

Fig 15. Profile of rates energy dissipation components at the minimum point of Ndf .

The instantaneous energy dissipation components of gasparticles drag forces dr and collisional interactions co at the minimum value point are shown in Fig. 17 at two computational cells. The mean energy dissipation components of drag forces and collisional interactions are 0.019 W/kg and 0.026 W/kg near the wall, and 1.06 W/kg and 0.419 W/kg at the center. Near the wall, the energy dissipation component of collisional interactions is larger than that of drag forces. It is reverse at the center, and the energy dissipation component of drag forces is larger than that of collisional interactions, indicating that the different controlling mechanisms of the formation of heterogeneous structures exist in the center and near the walls. If the energy dissipation component of gas-particles drag forces is smaller than the energy dissipation component of collisional interactions, the conditioned extreme value Eq. (24) simplifies to

εg , p, ug , us , θ

BEV BEV Ndf ≈ Ndf,CD = co

εd ∈R 1 , εc ∈R 2

Fig. 16. Instantaneous drag energy dissipation and collisional energy dissipation components at the minimum point of Ndf .

results in drag energy dissipation dr . The gas shear strain rates in the dense and dilute phases give rise to gas viscous energy dissipation vi . Fig. 15 shows the profile of instantaneous energy dissipation components at the minimum point. The simulated drag energy dissipation component dr and collisional energy dissipation component co are of the same order of magnitude. However, the gas viscous energy dissipation component vi is several orders of magnitude smaller than the drag and collisional energy dissipation components. Therefore, the authors assume that the gas viscous energy dissipation component vi is negligible in the determination of the minimum value point of energy dissipation rate Ndf of the computational cell. Since the order of magnitude of gas viscous energy dissipation component is much less than drag energy dissipation and collisional energy dissipation components, Eq. (17) can further be simplified to the following equation ε , p, u , u , θ BEV Ndf ≈ ( dr + co )|εg ∈R1 , εg c ∈Rs 2 = minimum d

(24)

The distributions of instantaneous energy dissipation components dr and co at the minimum value point are shown in Fig. 16 in the riser. The high slip velocity between the gas phase and the solid phase gives large energy dissipation component dr . The collisional energy dissipation component co is large at the high solid volume fractions. Generally, the energy dissipation component dr is large near the inlet because of high inlet gas velocities. Furthermore, it is high at the center due to large gas velocity. Near the walls the energy dissipation component co is large because of particles accumulation at the walls.

= minimum

dr < co

(25)

The minimization of energy dissipation rate is equivalent to the minimization of energy dissipation component of collisional interactions of heterogeneous structures. It appears that dr <co must be satisfied for Eq. (24) to be a valid approximation since the error in Eq. (25) is of the order 0(dr ). This indicates the interactions of collisions are dominant, and it is known as a collision-dominant (CD) heterogeneous structure. In the granular homogenous cooling system (HCS) with high solids volume fractions, the particles are agglomerated because the granular temperature decreases due to particle–particle collisions (Haff et al., 1983; McMillan et al., 2013; Yin et al., 2013). The formation of heterogeneous structure could arise when the collisions between the particles are sufficiently inelastic. When the energy dissipation component of gas-particles drag forces is larger than the energy dissipation component of collisional interactions, the energy dissipation component of hydrodynamic interactions is dominant. The conditioned extreme value Eq. (24) simplifies to

εg , p, ug , us , θ

BEV BEV Ndf ≈ Ndf,HD = dr 

εd ∈R 1 , εc ∈R 2

= minimum

dr > co (26)

It appears that dr >co is satisfied for Eq. (24) to be a valid approximation since the error in Eq. (26) is of the order 0(co ), and it is known as a hydrodynamic-dominant (HD) heterogeneous structure. The formation of clusters arises from the hydrodynamic forces of heterogeneous structures. Wylie & Koch (20 0 0) found that the formation of clusters was driven by the energy dissipation due to interactions surrounding gas and particles. The dilute clusters with low solids volume fraction increase slip velocity between the gas and particles, thereby, tend to be destroyed due to the increasing gas flow within the cluster (Moran and Glicksman, 2003). Interactions between gas and particles promote energy dissipation

14

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

Fig. 17. Simulated instantaneous drag energy dissipation and collisional energy dissipation components at the minimum point.

component by drag forces due to the velocity difference between the gas phase and the solid phase (Eaton and Fessler, 1994). From Eq. (16), Eq. (26) of the HD controlling mechanism is

 Ndf,HD =



εg , p, ug , us , θ nd FdUd + nc FcUc + ni FiUi (1 − f )  = minimum  εg ρg + (1 − εg )ρs

(27) Replacing solid mass flux Gs and superficial gas velocity Ug with p, ug and θ , Eq. (27) becomes

 Ndf,HD =

nd FdUd + nc FcUc + ni FiUi (1 − f ) ε g ρg + ( 1 − ε g ) ρs

εg , Ug , Gs   = minimum  (28)

It is known as the stability condition used in the EMMS drag models (Li and Kwauk, 2003; Nikopouls et al., 2010; Shah et al., 2011; Schneiderbauer et al., 2013; Zeneli et al., 2015). This indicates that the stability condition Eq. (28) used in the EMMS drag models is one of the controlling mechanisms of the DCSD drag model from Eq. (24) for flow of heterogeneous structures. The controlling mechanisms of CD and HD are marked as I and II in Fig. 17. When the energy dissipation component dr is the same order as the collisional energy dissipation component co , the heterogeneous structure is compromised by the gasparticles interactions and the interactions of collisions of particles, and identifies a collision-hydrodynamic-dominant (CHD) heterogeneous structure. The heterogeneous structures are formed with the variation of drag energy dissipation component and collisional energy dissipation component in the computational cell. In order to identify the contributions of the hydrodynamic interactions and interactions of collisions on heterogeneous structures, the ratio of energy dissipations dr /co defines the controlling mechanism index of heterogeneous structure.

ι =

⎧ ⎨ HD

dr /co > A1

CHD

A1 ≥ dr /co ≤ A2

CD

≥ dr /co ≥ A2

⎩

(29)

where the subscript i represents HD, CHD and CD, respectively. A1 and A2 are the base line values. Here, we assume that the base lines of A1 and A2 are 10.0 and 0.1, respectively. When the ratio dr /co is larger than 10.0, the energy dissipation component by drag forces is one order of magnitude larger than the collisional energy dissipation component, and the heterogeneous structure identifies a HD. When the ratio dr /co is less than 0.1, the existence of clusters is dominated by energy dissipation component co , representing a CD heterogeneous structure. When the ra-

Fig. 18. Distributions of controlling mechanism index i (i=CD, CHD, HD and HF) of clusters.

tio dr /co falls in between 0.1 and 10.0, the existence of clusters is compromised by the hydrodynamic interactions and the interactions of collisions of particles, representing a CHD heterogeneous structure. Thus, the controlling mechanism index i represents an occurrence of the controlling mechanism i in the computational cells. The distributions of instantaneous controlling mechanism index are shown in Fig. 18 in the riser. The variation of instantaneous controlling mechanism index suggests that the formation of clusters is dominated by hydrodynamics interactions and interactions of collisions of particles in the riser. The fraction of the controlling mechanism ψ i defines the ratio of the occurrence number of i to the total number N of computational cells of the riser.

ψι =



i

N

(30)

where the subscript i represents HD, CHD, CD and HF, respectively. The sum of the fractions of ψ i is 1.0. The instantaneous fractions of the controlling mechanism of heterogeneous structures are shown in Fig. 19 in the riser. The mean fractions of the controlling mechanism of CD, HD and CHD are 0.038, 0.115 and 0.298 of the riser, suggesting that the controlling mechanism CD of heterogeneous structures is lowest, and the next in importance is the controlling mechanism HD. The fractions of the controlling mechanism CHD of heterogeneous structures is largest, indicating the heterogeneous structures arise as a result of two controlling mechanisms operat-

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

15

Fig. 21. Distribution of solid volume fractions of dense phase. Fig. 19. Instantaneous fraction of controlling mechanism components of the riser.

Fig. 20. Profile of predicted cluster diameters with solid volume fractions. Fig. 22. Profile of granular temperatures of clusters and dispersed particles.

ing in parallel for the inelastic collisions of particles and the hydrodynamic interactions. 4.6. Predicted cluster properties from DCSD drag model From simulated instantaneous cluster diameters, the timeaveraged cluster diameters are shown in Fig. 20 with the change of solid volume fractions. The large variation of cluster diameters originates from the change of computational cell parameters. Generally, the cluster diameters are small at the low solid volume fractions. The empirical correlations of cluster diameters proposed by Gu and Chen (1998) and Harris et al. (2002) give a monotonous distribution with the change of solid volume fractions. On the other hand, the predictions using the CD-Lab model (Schneiderbauer et al., 2013) increase, reach maximum and then decrease with the increase of solid volume fractions. Both the DCSD drag model and the CD-Lab model show similar trends of predicted cluster diameters. The difference between the simulations and calculations using empirical correlations is obvious. Several reasons lead to the difference between simulations and experiments. Firstly, the DCSD model assumes one-dimensional flow of clusters. A three-dimensional model is needed to describe flow of 3-D clusters. Another factor is that the measured mean vertical length was taken as the size of clusters because the shape of clusters varies considerably (Harris et al., 2002). Further the difference is introduced by experimental measurements because the measured clusters are at or near the wall of a riser. Studies of cluster size in the riser core are complexity (Harris et al., 2002). Hence, the complete understanding of distribution of diameter of clus-

ters requires measurements and development of 3-D drag model of heterogeneous structures. The distribution of time-averaged cluster solid volume fraction is shown in Fig. 21 with the change of solid volume fractions in the computational cells. The calculated cluster solid volume fractions from the empirical correlations proposed by Gu and Chen (1998) and Harris et al. (2002) are also plotted with solid volume fractions. The cluster existence was identified by Soong et al. (1993) to detect cluster solid volume fraction. The criteria were the solid volume fraction in a cluster must be two times the standard deviation greater than the time averaged solid fraction at that measuring position (Sharma et al., 20 0 0). A similar definition of the identification of clusters was proposed by Lints and Glicksman (1993) based upon experimental measurements. These definitions identify cluster solid volume fraction based upon a threshold criterion for solid volume fraction, indicating the empirical correlation of cluster solid volume fraction relates to the threshold from experimental measurements. The distributions of time-averaged granular temperatures of dense phase and dilute phase are shown in Fig. 22 as a function of load ratio. Generally, the granular temperature of dense phase is lower than that of dilute phase. From experimental measurements using a laser Doppler velocimetry, Breault et al. (2005) obtained granular temperatures of clusters and dispersed particles in a fluidized bed. They found that the granular temperatures of cluster are different from the dispersed particles granular temperature. All downward moving clusters have a mean granular temperature of 0.032 (m/s)2 , while for the upward moving clusters it has a mean

16

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126

granular temperature of 0.12 (m/s)2 . The granular temperature of dispersed particles is from 0.4 to 0.9 (m/s)2 , indicating the granular temperature of dispersed particles is an order of magnitude larger than that for clusters. The simulated mean granular temperatures of clusters and dispersed particles are 0.00229 (m/s)2 and 0.0666 (m/s)2 , respectively. The measured granular temperature is an order of magnitude larger than that simulated by DCSD because there are limitations in the DCSD model with its one dimensionality. This is a good basis for the start of the simulations of heterogeneous structures but it certainly is an approximation, indicating the DCSD drag model needs improvement further. 5. Conclusions A DCSD drag model is proposed using the minimization of energy dissipation rate to take the hydrodynamic interactions and collisional interactions into account in gas-solids CFB risers. The heterogeneous structure parameters are predicted using the conditioned extreme value equation which distinguishes from heterogeneous flow and homogeneous flow in computational cells. The DCSD drag coefficient is predicted by means of the set of transient nonlinear equations with the BEV theory as a function of ten heterogeneous structure parameters according to five computational cell parameters. The simulated low intermittency factor indicates that the occurrence of heterogeneous structures is intermittent in the computational cells. Three controlling mechanisms named as CD, HD and CHD mechanisms are defined according to the energy dissipation component of hydrodynamic interactions and the collisional interaction. The formation of clusters results from the compromise in competition of hydrodynamic interactions and interactions of collisions of particles by CHD mechanism in the riser. The comparison between the numerical simulations and the calculations using empirical correlations suggests that the DCSD drag model needs worth further studying in the future. In view of flow of heterogeneous structures, the general expressions for the conservation equations of momentum and granular temperature of particles in the dilute phase and clusters in the dense phase are required. However, we still lack an understanding of the relation between the velocity fluctuations of the gas and the particles in the dense and dilute phases. More research is needed to render for an improvement of the formulations of dispersed particles in the dilute phase and clusters in the dense phase. Acknowledgment This work is supported by the Natural Science Foundation of China through grant nos. 91752115 and 51776059. Reference Agrawal, K., Loezos, P.N., Syamlal, M., Sundaresan, S., 2001. The role of meso-scale structures in rapid gas-solid flows. J. Fluid Mech. 445, 151–185. ANSYS, Inc., Ansys fluent theory guide, 2011. Arastoopour, H., Gidaspow, D., Abbasi, E., 2016. Computational Transport Phenomena of Fluid-Particle System. Springer, New York. Basu, P., 2015. Circulating Fluidized Bed boilers: Design, Operation and Maintenance. Springer, Cham. Breault, R.W., 2006. A review of gas–solid dispersion and mass transfer coefficient correlations in circulating fluidized beds. Powder Technol. 163, 9–17. Breault, R.W., Ludlow, C.J., Yue, P.C., 2005. Cluster particle number and granular temperature for cork particles at the wall in the riser of a cfb. Powder Technol. 149, 68–77. Brereton, C.M.H., Grace, J.R., 1993. Micro-structural aspects of the behavior of circulating fluidized beds. Chem. Eng. Sci. 48, 2565–2572. Cahyadi, A., Anantharaman, A., Yang, S., Karri, S.B.R., Findlay, J.G., Cocco, A.R.O, Chew, J.W., 2017. Review of cluster characteristics in circulating fluidized bed (CFB) risers. Chem. Eng. Sci. 158, 70–95. Cloete, J.H., Cloete, S., Municchi, F., Radl, S., 2018. Development and verification of anisotropic drag closures for filtered two fluid models. Chem. Eng. Sci. 192, 930–954.

Cocco, R., Shaffer, F., Hays, S.B.R., Karri, T., Knowlton, R., 2010. Particle clusters in and above fluidized beds. Powder Technol. 203, 3–11. Eaton, J.K., Fessler, J.R., 1994. Preferential concentration of particles by turbulence. Int. J. Multiph. Flow 20, 169–209. Fox, R.O., 2014. On multiphase turbulence models for collisional fluid-particles flows. J. Fluid Mech. 742, 368–424. Gao, X., Tingwen, L., Rogers, W.A., 2018. Assessment of mesoscale solid stress in coarse-grid tfm simulation of geldart a particles in all fluidization regimes. AIChE J. 64, 3565–3581. Gidaspow, D., 1994. Multiphase Flow and fluidization: Continuum and Kinetic Theory Descriptions. Academic Press, New York. Gidaspow, D., Jiradilok, V., 2009. Computational techniques: Energy science, Engineering and Technology Series. NOVA, New York. Goldhirsch, I., Zanetti, G., 1993. Clustering instability in dissipative gases. Phys. Rev. Lett. 70, 1619–1622. Grace, J.R., Avidan, A.A., Knowlton, T.M., 1997. Circulating Fluidized beds, Blackie Academic & Professional. Springer. Grace, J.R., Tuot, J., 1979. A theory for cluster formation in vertically conveyed suspensions of intermediate density. Trans. Ind. Chem. Eng. 57, 49–54. Gu, W.K., Chen, J.C., 1998. A model for solid concentration in circulating fluidized beds. In: Fan, L.S., Knowlton, T.M. (Eds.), Proceedings of the 9th Engineering Foundation Conference on Fluidization, Durango, CO. Engineering Foundation, pp. 501–508 NY. Gustavsson, K., Vajedi, S., Mehlig, B., 2014. Clustering of particles falling in a turbulent flow. Phys. Rev. Lett. 112, 214501. Haff, P.K., 1983. Grain flow as a fluid mechanical phenomenon. J. Fluid Mech. 134, 401–948. Harris, A.T., Davidson, J.F., Thorpe, R.B., 2002. The prediction of particle cluster properties in the near wall region of a vertical riser. Powder Technol. 127, 128–143. Huilin, L., 2017. Kinetic Theory For Dense Fluid-Solid Flow. Science Press. Igci, Y., Andrews, A.T., Sundaresan, S., Pannala, S., O’Brien, T., 2008. Filtered two-fluid models for fluidized gas-particle suspensions. AIChE J. 54, 1431–1448. Johnson, P.C., Jackson, R., 1987. Frictional–collisional constitutive relations for granular materials, with application to plane shearing. J. Fluid Mech. 176, 67–93. Kotz, S., Nadarajah, S., 20 0 0. Extreme Value distributions, Theory and Applications. Imperial College Press, London. Li, D., Wang, S., Liu, G., Huilin, L., Jiang, X., Ming, T., Zhenjie, L., 2018. A dynamic cluster structure-dependent drag coefficient model applied to gas-solid risers. Powder Technol. 325, 381–395. Li, J., Kwauk, M., 2003. Exploring complex systems in chemical engineering: the multi-scale methodology. Chem. Eng. Sci. 58, 521–535. Li, T., Pougatch, K., Salcudean, M., Grecov, D., 2008. Numerical simulation of horizontal jet penetration in a three-dimensional fluidized bed. Powder Technol. 184, 89–99. Lints, M., Glicksman, L.R., 1993. The structure of particle clusters near the wall of a circulating fluidized bed. AIChE Symp. Ser. 89, 35–47. Luding, S., Huthmann, M., McNamara, S., Zippelius, A., 1998. Homogeneous cooling of rough dissipative particles: theory and simulations. Phys. Rev. E 58, 3416–3425. Manyele, S.V., Parssinen, J.H., Zhu, J.X., 2002. Characterizing particle aggregates in a high-density and high-flux cfb riser. Chem. Eng. J. 88, 151–161. McKeen, T., Pugsley, T., 2003. Simulation and experimental validation of a freely bubbling bed of fcc catalyst. Powder Technol. 129, 139–152. McMillan, J., Shaffer, F., Gopalan, B., Chew, J.W., Hrenya, C.M., Hays, R., Karri, S.B.R., Cocco, R., 2013. Particle cluster dynamics during fluidization. Chem. Eng. Sci. 100, 39–51. Molerus, O., Wirth, K.E., 1997. Heat Transfer in Fluidized Beds. Springer, Dordrecht. Mouallem, J., Niaki, S.R.A., Norman, C.C., Miliolo, C.C., Milioli, F.E., 2019. Macro-scale effects over filtered and residual stresses in gas-solid riser flows. Chem. Eng. Sci. 195, 553–564. Moran, J.C., Glicksman, L.R., 2003. Experimental and numerical studies on the gas flow surrounding a single cluster applied to a circulating fluidized bed. Chem. Eng. Sci. 58, 1879–1886. Nikolopoulos, A., Papafotiou, D., Nikolopoulos, N., Grammelis, P., Kakaras, E., 2010. An advanced emms scheme for the prediction of drag coefficient under a 1.2MWth cfbc isothermal flow-Part I: numerical formulation. Chem. Eng. Sci. 65, 4080–4088. Ozarkar, S.S., Yan, X., Wang, S., Milioli, C.C., Milioli, F.E., Sundaresan, S., 2015. Validation of filtered two-fluid models for gas–particle flows against experimental data from bubbling fluidized bed. Powder Technol. 284, 159–169. Parmentier, J.F., Simonin, O., Delsart, O., 2012. A functional subgrid drift velocity model for filtered drag prediction in dense fluidized bed. AIChE J. 58, 1084–1098. Sarkar, A., Miloioli, F.E., Ozarkar, S., Li, T., Sun, X., Sundaresan, S., 2016. Filtered sub-grid constitutive models for fluidized gas-particle flows constructed from 3-D simulations. Chem. Eng. Sci. 152, 443–456. Schneiderbauer, S., Puttinger, S., Pirker, S., 2013. Comparative analysis of subgrid drag modifications for dense gas-particle flows in bubbling fluidized beds. AIChE J. 59, 4077–4099. Schneiderbauer, S., Pirker, S., 2014. Filtered and heterogeneity-based subgrid modifications for gas–solid drag and solid stresses in bubbling fluidized beds. AIChE J. 60, 839–854. Shah, M.T., Utikar, R.P., Tade, M.O., Pareek, V.K., 2011. Simulation of gas-solid flows in riser using energy minimization multiscale model: effect of cluster diameter correlation. Chem. Eng. Sci. 14, 3291–3300.

J. Xiaoxue, L. Dan and W. Shuyan et al. / International Journal of Multiphase Flow 122 (2020) 103126 Sharma, A., Tuzla, K., Matsen, J., Chen, J., 20 0 0. Parametric effects of particle size and gas velocity on cluster characteristics in fast fluidized beds. Powder Technol. 111, 114–122. Shuai, W., Guangbo, Z., Guodong, L., Huilin, L., Feixiang, Z., Tianyu, Z., 2014. Hydrodynamics of gas–solid risers using cluster structure-dependent drag model. Powder Technol. 254, 214–227. Subbarao, D., 2010. A model for cluster size in risers. Powder Technol. 199, 48–54. Soong, C.H., Tuzla, K., Chen, J.C., 1993. Identification of particle clusters in circulating fluidized beds. In: Avidan, A.A. (Ed.), Circulating Fluidized Bed Technology IV. Engineering Foundation, New York, pp. 615–620. Syamlal, M., O’Brien, T.J., 2003. Fluid dynamic simulation of O3 decomposition in a bubbling fluidized bed. AIChE J. 49, 2793–2801. Syamlal, M., Rogers, W., O’Brien, T.J., 1993. MFIX documentation theory guide. U.S. Department of Energy, Office of Fossil Energy, Morgantown, West Virginia Technical note. Technical report. Wei, F., Yang, G., Jin, Y., Yu, Z., 1995. The characteristics of clusters in a high density circulating fluidized bed. Can. J. Chem. Eng. 73, 650–655. Wei, F., Lin, H., Cheng, Y., Wang, Z., Jin, Y., 1998. Profiles of particle velocity and solids fraction in a high-density riser. Powder Technol. 100, 183–189. Wen, C.Y., Yu, Y.H., 1966. Mechanics of fluidization. Chem. Eng. Progr. Symp. Ser. 62, 100–111.

17

Wenming, L., Hongzhong, L., Qingshan, Z., 2017. Modeling the hydrodynamics of downer reactor based on the meso-scale structure. Powder Technol. 314, 367–376. Wu, R.L., Lim, C.J., Grace, J.R., Brereton, C.M.H., 1991. Instantaneous local heat transfer and hydrodynamics in a circulating fluidized bed. Int. J. Heat Mass Transf. 34, 2019–2027. Wylie, J.J., Koch, D.L., 20 0 0. Particle clustering due to hydrodynamic interactions. Physics of Fluids 12, 964–970. Yates, J.G., Lettieri, P., 2016. Fluidized-bed reactors: Processes and Operating Conditions. Berlin: Springer. Yerushalmi, J., Cankurt, N.T., Geldart, D., Liss, B., 1976. Flow regimes in vertical gas— solid contact systems. AIChE Symp. Ser. 74, 1–13. Yin, X., Zenk, J.R., Mitran, P.P., Hrenya, C.M., 2013. Impact of collisional versus viscous dissipation on flow instabilities in gas-solid systems. J Fluid Mech. 727, 1–12. Zeneli, M., Nikolopoulos, A., Nikolopoulos, N., Grammelis, P., Kakaras, E., 2015. Application of advanced coupled emms-tfm model to a pilot scale cfb carbonator. Chem. Eng. Sci. 138, 482–498. Zhang, D.Z., VanderHeyden, W.B., 2002. The effects of mesoscale structures on the macroscopic momentum equations for two-phase flows. Int. J. Multiph. Flow 28, 805–822.