Comparison of different drag models in CFD-DEM simulations of spouted beds

PTEC-14825; No of Pages 18 Powder Technology xxx (2019) xxx

Contents lists available at ScienceDirect

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

Comparison of different drag models in CFD-DEM simulations of spouted beds Filippo Marchelli a,b, Qinfu Hou b,⁎, Barbara Bosio c, Elisabetta Arato c, Aibing Yu b,d a

Faculty of Science and Technology, Free University of Bozen-Bolzano, 39100 Bolzano, Italy ARC Research Hub for Computational Particle Technology, Department of Chemical Engineering, Monash University, Clayton, VIC 3800, Australia Process Engineering Research Team, Department of Civil, Chemical and Environmental Engineering, University of Genova, 16145 Genova, Italy d Center for Simulation and Modelling of Particulate Systems, Southeast University - Monash University Joint Research Institute, Suzhou 215123, PR China b c

a r t i c l e

i n f o

Article history: Received 10 July 2019 Received in revised form 8 October 2019 Accepted 12 October 2019 Available online xxxx Keywords: Eulerian-Lagrangian approach Fluidisation Gas-solid exchange coefficient Spouted bed User-defined function

a b s t r a c t Spouted beds are commonly simulated through the Computational Fluid Dynamics – Discrete Element Method approach. The choice of the drag model is still a matter of debate, as they feature peculiar operative conditions. In this work, we simulated two spouted beds containing Geldart-D particles. We tested seven drag models: three are classic models, while four are developed through advanced computational techniques. The results indicate that the key variable is the ratio between the operative and the minimum spouting gas velocity (u/ums). At u = ums only the Gidaspow model can always predict fluidisation, but at low u/ums values the Beetstra model is the best compromise. For higher values, the Rong and Di Felice models behave better, while the others overestimate the particles' velocity. These results can be useful to identify the best performing model and show there is a need for more appropriate models for spouted beds. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Spouted beds (SBs) are a type of gas-solid contactor. They share some features with traditional fluidised beds, but differ in the feeding of the fluidising agent, which, in the classic configuration, only happens through a single nozzle in the centre of the base. Thus, the particles undergo trajectories that show three distinct zones: the spout, the fountain and the annulus. Fig. 1, taken from a simulation of a SB, provides a visual depiction of these three zones. The thus-originated vigorous dynamic flow promotes high mass and energy transfer rates and makes SBs adequate to process coarse and irregular particles. As a matter of fact, SBs were initially proposed in Canada halfway through last century [1] with the aim of drying wheat particles. Throughout the years, researchers exploited SBs' advantages to employ them in several other applications, including: drying of seeds or other food products [2,3], coating [4], mixing [5], polymerisation [6], combustion [7], pyrolysis [8,9] and gasification [10–12] of coal and biomass, torrefaction [13] and electrochemical reactions [14]. For some of these applications, alternative configurations have also been proposed, such as spout-fluidised beds [15], multiple-stage spouted beds [6,9], spouted beds with draft tubes [16] and fountain confiners [17], and prismatic spouted beds [18]. Some of these have been employed at the industrial scale [5,19];

⁎ Corresponding author. E-mail address: [email protected] (Q. Hou).

the German company Glatt currently produces prismatic units for various applications [20]. Numerical studies of SBs are also quite abundant. Much effort has particularly been devoted to the simulation of the fluid dynamics of SBs, which still pose some modelling problems. Both the two-fluid model (TFM) and the coupled computational fluid dynamics – discrete element method (CFD-DEM) have been applied for this purpose. An established simulation methodology could provide reliable scale-up criteria, the lack of which has so far hindered an extensive industrial application of SBs. Several authors have recently reviewed the efforts done in the numerical modelling of SBs [21–24]. One major point is the lack of agreement between researchers on the most suitable drag model. Drag is the main cause of momentum exchange between a fluid phase and solid particles in contact with each other. The drag force on an isolated particle can be easily calculated. Conversely, the drag force acting in dense systems has a complex dependency on the porosity and on the particle Reynolds number. Finding an equation for this force has attracted the efforts of researchers for several decades now. At first, the studies were mainly based on experimental observations, with the force empirically estimated from the bed pressure drop. In this framework, the most popular models are those by Ergun [25], Wen-Yu [26] and Di Felice [27], which are still considered reasonably reliable and are widely applied. More recently, advanced computational techniques, such as Lattice-Boltzmann simulations, enabled a more direct estimation of the drag force. Several authors thus proposed new drag models, more complex than the earlier ones.

https://doi.org/10.1016/j.powtec.2019.10.058 0032-5910/© 2019 Elsevier B.V. All rights reserved.

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

2

F. Marchelli et al. / Powder Technology xxx (2019) xxx

Fig. 1. Spouting regime in a pseudo-2D spouted bed (particles coloured by the solid volume fraction).

Numerical simulations have been the tool for researchers to compare the performance of different drag models [28–35]. However, most of the studies are applied to traditional fluidised beds. Compared to them, SBs commonly feature coarser particles, higher particle Reynolds numbers, and distinct zones with quite constant porosities, so the same conclusions cannot be drawn a priori. As reported by Moliner et al. [23], the Gidaspow drag model is by far the most employed for TFM and CFD-DEM simulations of SBs, even though the choice is often not explained. Comparison of different drag models in SBs are however scarce; Table 1 summarises their findings. From the Table, it is evident that there is no agreement regarding the most suitable drag model for SBs. The only model that all authors tested is the one developed by Gidaspow, and some of them deemed it as better performing than others. It is a combination of the Wen-Yu and Ergun models, with a switch based on the porosity. In these early models, the drag force has quite simple equations, and mainly depends on the porosity and particle Reynolds number, with constant exponents. This was partially corrected by Di Felice, who showed that the exponent of the porosity depends on the particle Reynolds number. More recently, researchers took advantage of Lattice-Boltzmann simulations to propose more complex and accurate equations. For example, Rong and colleagues [59] showed that the Di Felice model can be updated by adding a dependence on the porosity in the exponent, while others proposed entirely new and more complex formulations. Indeed, as Table 1 shows, when these models are included in the comparison, they perform better than the classic ones. The superiority of one of these newer models over another is however unclear. Moreover, their application for SBs has so far been limited and needs to be further investigated. This will allow understanding whether the fluid dynamic

conditions of SBs fall within the range of applications of these drag models. The Discrete Element Method was developed by Cundall and Strack [60] and first coupled with CFD by Tsuji and colleagues [61,62]. The CFDDEM approach is widely applied for SBs, and is reportedly more precise than the TFM [63–65]. In our previous work [57], we used CFD-DEM to perform a sensitivity analysis on a SB. The analysis highlighted the importance of including a turbulence and a Magnus lift model, and the severe impact of the drag models on the particles' behaviour. In this work, the focus is put on the drag force and more models are tested, including models developed by means of Lattice Boltzmann or direct numeric simulations. In summary, seven models are included in the study: three (Wen-Yu, Gidaspow, Di Felice) are classic choices for this kind of application, two (Beetstra, Koch-Hill) have recently deemed as good choices for SBs, and two (Rong, Tenneti) have recently been applied to fluidised beds with good results, but not yet to SBs. Moreover, two different SBs are simulated to study a wider range of operative conditions and possibly widen the validity of our findings. The first device [66,67] is the SB already employed in the previous work: it is a pseudo-2D unit, as its depth is much lower than its width and height. The other is a cylindrical SB developed by Olazar and colleagues [68,69], with identical width and depth, and hence fully 3D. Both reactors contain monosized coarse Geldart [70] particles. To the knowledge of the authors, this is the first work comparing the performance of the same drag models in two different SBs, and it is the first time the Rong and Tenneti drag models are applied to SBs. 2. Governing equations Nowadays, the CFD-DEM approach is a popular tool among researchers. Its core principle is the coupled resolution of the continuum balance equations for the gas phase, and of the discrete equations of motion of each solid particle for the solid phase. In this work, we applied the CFD-DEM approach based on its formulation in the commercial code Ansys Fluent 19.2, using the Dense Discrete Phase Model (DDPM), which is the setting in Fluent that allows including a discrete phase. For the sake of brevity, only the main equations are reported, and the reader is referred to the cited Fluent guides for more information. 2.1. Gas phase The behaviour of the gas phase q is predicted with the local mass and momentum balance equations [71–73], which Table 2 contains. ! In Eq. 2, R pq is the momentum exchange between the gas (q) and solid (p) phases:
 ! ! ! R pq ¼ K pq u p − u q

ð8Þ

in which Kpq is the gas-solid exchange coefficient, which depends on the

Table 1 Works from the literature comparing drag models in spouted beds. Geometry

Modelling approach

Tested models

u/ums

Ref.

Cylindrical spouted bed Cylindrical spouted bed Cylindrical spouted bed Conical spouted bed Prismatic spouted bed Flat-bottomed spout-fluid bed Pseudo-2D flat-bottomed spouted bed Pseudo-2D pyramidal spouted bed Pseudo-2D pyramidal spouted bed Pseudo-2D pyramidal spouted bed

TFM TFM TFM TFM CFD-DEM CFD-DEM CFD-DEM CFD-DEM CFD-DEM TFM

Richardson and Zaki [36], Gidaspow [37], Syamlal-O'Brien [38], Di Felice [27], Arastoopour [39] Gidaspow, Gao et al. [41], Zhang and Reese [42] and three other self-developed hybrid models Gidaspow and Du Plessis [44] Self-developed model and Huilin-Gidaspow model [46] Gidaspow, Di Felice, Koch and Hill [48], Beetstra [49] Gidaspow, Du Plessis and Masliyah [51], Van Der Hoef [52] and Beetstra Wen-Yu, Gidaspow, Syamlal-O'Brien Dahl and Hrenya [55], Syamlal-O'Brien, Ding-Gidaspow, Lathouwers-Bellan, Koch-Hill, Beetstra Wen-Yu, Gidaspow, Syamlal-O'Brien Syamal-O'Brien, Di Felice, Wen-Yu, Gibilaro, Gidaspow and Huiling - Gidaspow

1.1, 1.2, 1.3. Not stated. 1.0, 1.1, 1.2, 1.3. Not stated. Not stated. Not stated. Not stated. 1.74 1.74 1.74

[40] [43] [45] [47] [50] [53] [54] [56] [57] [58]

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

F. Marchelli et al. / Powder Technology xxx (2019) xxx

3

Table 2 Equations solved for the gas phase. Name

Equation

Continuity equation

∂ ! ðaq ρq Þ þ ∇∙ðaq ρq u q Þ ¼ 0 ∂t ! ∂ ! ! ! ! ! ! ðaq ρq u q Þ þ ∇∙ðaq ρq u q u q Þ ¼ −aq ∇p þ ∇∙t q þ aq ρq g þ R pq ∂t ! ! ! ! ! !T 2 ! t q ¼ − ðρq kq þ ρq νt;q ∇∙ U q Þ I þ ρq νt;q ð∇ U q þ ∇ U q Þ 3 ∂ ! ðaq ρq kq Þ þ ∇∙ðaq ρq u q kq Þ ¼ ∇∙½aq ρq ðυq þ υT;q Þ∇kq  þ aq Gk;q −aq ρq εq ∂t εq ∂ ! ðaq ρq εq Þ þ ∇∙ðaq ρq u q εq Þ ¼ ∇∙½aq ρq ðυq þ υT;q Þ∇εq  þ aq ð1:44Gk;q −1:92ρq εq Þ kq ∂t ∂u j Gk;q ¼ −ρq uq;i uq; j ∂xi

Momentum balance equation Stress-strain tensor Turbulence kinetic energy Turbulence dissipation rate Production of turbulence kinetic energy Turbulent viscosity

kq εq

drag force (discussed in detail in the following section):
 ! ! ! K pq u p − u q ¼ ag as ρp F d

ð9Þ

under the assumption that, in each cell: np mp ¼ as ρp

ð10Þ

In the previous equations, ui is the velocity of phase i, ai is the volume fraction of the phase i, ρi is the density of phase i, Fd is the drag force, np is the number of particles per cell and mp is the mass of a particle. In CFD simulations, a correct description of the gas phase behaviour often relies on the inclusion of a turbulence model. As a matter of fact, in our previous work [57] we showed that neglecting turbulence modelling in a SB resulted in an overestimation of the particle velocity. Conversely, the common two-equation turbulence models produce very similar results. Other more complex approaches exist, such as the Large Eddy Simulation (LES) method, which can be coupled with the DEM [74]. However, this coupling is currently not feasible in Fluent. Furthermore, CFD-DEM simulations of coarse particles must employ coarse grids, which may limit the accuracy of turbulence models, especially near walls. In the present work, turbulence is considered through the inclusion of the standard k-ɛ model, in its dispersed form. The equations of the model can be found in Table 2. 2.2. Solid phase 2.2.1. Equations of motion The DEM calculates the trajectory of each of the particles that the solid phase comprises. To do so, the Newtonian equations of motion must be solved for each particle. These are: mp

Ip

! ! ! ! ! dup !ρp −ρ ! ¼ mp g þ F d þ F vm þ F pg þ F ml þ F c ρp dt

! dω p ρ
mp 5 ! ! ¼ C ω  Ω ∙Ω 2 2 dt

(2) (3) (4) (5) (6) (7)

2

υT;q ¼ 0:09

(1)

ð11Þ

ð12Þ

! ! ! ! ! In the first equation, the terms F d ; F vm ; F pg ; F ml and F c represent respectively the drag, virtual mass, pressure gradient, Magnus lift and contact forces. Being the focus of this work, the description of the drag force is treated in detail in the following Section. The meaning and calculation of the other forces and parameter are instead summarised in Table 3. The contact force is the resulting force from all the friction and collision forces acting on a particle at a given time. Here, collisions are modelled through a linear spring-dashpot model, which has produced good results in previous works [57] and is common in the literature about CFD-DEM modelling of SBs [15,54,66,75]. It requires the

specification of two parameters only: the spring constant (K) and the restitution coefficient for the dashpot term (η). It also has the advantage that the produced results are insensitive to the value of K, provided it is higher than 102–103 N/m. The other common approach to compute collision is the Hertzian one [50,53], which has a high computational cost. We wish to underline that the Magnus lift force, in the formulation by Rubinow and Keller [76], is included since our previous work showed that it plays a major role in SBs. This is coherent with recent findings of other authors, who showed that the Magnus lift force can become even larger than the drag force for high particle Reynolds number and porosities [77,78]. The equation of the pressure gradient force (Eq. 14) is different from the standard one and is based on the Fluent code. It relies on the hypothesis that diffusive and source terms are negligible, and so the pressure gradient can be replaced with the velocity gradient. The pressure gradient and virtual mass force have a very low magnitude in these simulations, as the density of the fluid phase is three orders of magnitude lower than that of the particles. 2.2.2. Drag force Drag is the main cause of momentum exchange between the gas and solid phases. As already explained in the Introduction, its modelling has a long history and is still not a completely established practice. To assess the behaviour of different drag models in CFD-DEM simulations of SBs, we tested several of them. In Fluent, the drag force has this general form: μ g
! !  ! F d ¼ mp u−up D 2 ρp dp

ð31Þ

where dp is the particle diameter, μg is the dynamic viscosity of the gas and D is the drag function. The dimensionless term D differs for each drag model, but it is always a function of the porosity and of the particle Reynolds number. In these simulations, the porosity coincides with the cell-based gas volume fraction ag. The particle Reynolds number is thus defined [73]: Rep ¼

! ! ρg dp  u p − u  μg

ð32Þ

In the formulations of several drag models, the particle Reynolds number is also multiplied by the porosity, so additional care must be taken upon implementing their equations in Fluent. The only drag models suitable for dense gas-solid systems that are available in Fluent are the Wen-Yu, Gidaspow and Syamlal-O'Brien models. To test other models, their equations must be supplied through user-defined-functions (UDF) written in the C language. The UDF provides a value for the function D. For a better comparison, in this work all the tested drag models were supplied through UFDs. The followed

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

4

F. Marchelli et al. / Powder Technology xxx (2019) xxx

Table 3 Closure equations for the solid phase. Name Virtual mass force Pressure gradient force Magnus lift force Rotational lift coefficient [76] Collision force for two contacting particles Equations of the spring-dashpot model

Friction force Rolling friction force Moment of inertia Relative particle-fluid angular velocity Rotational drag coefficient Rotational particle Reynolds number

Equation ! ! ! ρ ! ! dup F vm ¼ 0:5mp u p∇ u − ρp dt ! ρ ! ! F pg ¼ mp u p ∇ u ρp ! ! ! j u p− u j  ! ! 1 ! ! F ml ¼ AP C RL ρg ! ! ð u p − u Þ  ðω p −ω Þ 2 jω p −ω j ! jω p jdp C RL ¼ 2 ! ! 2j u p − u j  ! ! ! ! ! F 1 ¼ − F 2 ¼ Kδ þ γð u 12 ∙ e 12 Þ e 12 ! ! x 2− x 1 ! e 12 ¼ ! ! k x 2− x 1k ! ! δ ¼ k x 2 − x 1 k−ðr 1 þ r 2 Þ ! ! ! u 12 ¼ u 2 − u 1 m12 ln η γ ¼ −2 t coll rffiffiffiffiffiffiffiffiffi m12 t coll ¼ f loss K m1 m2 m12 ¼ m1 þ m2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 f loss ¼ π 2 þ ln η ! ! F friction ¼ μ f F normal ! ! F rolling ¼ μ rolling F normal π 5 ρ d Ip ¼ 60 p p ! 1 ! ! Ω ¼ ∇  u −ω p 2 6:45 32:1 C ω ¼ pffiffiffiffiffiffiffiffiffiffi þ Reω Reω ! ! 2 1 2 ∇  u −ω p Þdp Reω ¼ ρð 4μ

(13) (14) (15) (16) (17) (18) (19) (20) (21) (22) (23) (24) (25) (26) (27) (28) (29) (30)

procedure and the codes of the models are provided in the Supporting Information. Please note that the formulations of the following drag functions are only valid for spherical particles, and that the simulations only featured spherical particles. Recently, researchers also worked on developing drag models for non-spherical particles, which are obviously more complicated in their formulation and underlying principles [79–82]. While researchers have managed to perform CFD-DEM simulations of SBs featuring non-spherical particles [83,84], we decided to focus on the more fundamental case of spherical particles, for which it is more urgent to establish an agreement.

beds. Its formulation yields:

• Wen-Yu drag model

where CD is the drag coefficient proposed by Dallavalle:


8  0:687  −3:65 > 1 þ 0:15 ag Rep < 18 ag D¼ Rep asolid > : 150 2 þ 1:75 ag ag

for asolid N 0:2

ð34Þ

• Di Felice drag model Di Felice proposed the following drag model in 1994 [27]: D ¼ 0:75 C D Rep ag 2−β

The Wen-Yu [26] drag model results in the following D: CD ¼
 0:687  D ¼ 18 a−3:65 1 þ 0:15 ag Rep g

for asolid ≤ 0:2

ð35Þ

!2 4:8 0:63 þ pffiffiffiffiffiffiffiffiffiffiffiffi ag Rep

ð36Þ

ð33Þ

It is one of the first developed drag models, and it hypothesises that the coefficient of the porosity is constant at −3.65. This is a valid assumption for low and high Reynolds numbers, but it was later proved that there is an intermediate range in which the coefficient varies [27]. The model is considered mostly valid for dilute systems and is still often chosen due to its simple formulation. • Gidaspow drag model The Gidaspow drag model [85] is a combination of the Wen-Yu and the Ergun models, switching from the former to the latter when the porosity is lower than 0.8. This allows it to work well in both dense and in dilute regions. It is the most popular model for fluidised and spouted

Di Felice considered data from several systems and showed that the coefficient of ag should not be constant but should instead be a function of the particle Reynolds number, displaying a minimum at about Rep = 50. This is considered through β, which is a function of only Rep:  2 ! 1:5− log 10 ag Rep ¼ 3:7−0:65 exp − 2 

βDi Felice

ð37Þ

• Rong et al. model The model of Rong et al. [59] can be considered an enhancement of the Di Felice drag model. As a matter of fact, it uses the same D function, and only varies in the definition of β. Through several Lattice-Boltzmann simulations, Rong and colleagues showed that β cannot be considered a function of the sole particle Reynolds number, because the porosity

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

F. Marchelli et al. / Powder Technology xxx (2019) xxx

5

Table 5 List of reference parameters employed in the simulations (the contact parameters are the same for particle-particle and particle-wall interactions). Parameter

Pseudo-2D

Spring constant Restitution coefficient Friction coefficient Rolling friction coefficient Particle density Particle diameter Initial bed height Total number of particles Shape of particles

1000 N/m 0.9 0.3 0.03 2380 kg/m3 2.033 mm, monodispersed 7.5 or 10 cm 9399 or 15,590 Spherical

Cylindrical

2420 kg/m3 4 mm, monodispersed 20 cm 49,144

Beetstra, van der Hoef and Kuipers [49] also employed LatticeBoltzmann simulations of fluid flows through arrays of spheres. They tested several porosity values, and Reynolds numbers up to 1000. The model they thus developed leads to the following equation for D: D ¼ 18ag F B Fig. 2. Front and bottom view of the meshes of the pseudo-2D (a) and of the cylindrical (b) spouted beds.

where

plays a role as well: βRong

  2 !     1:5− log10 ag Rep ¼ 2:65 ag þ 1 − 5:3−3:5ag a2g exp − 2 ð38Þ

The Rong and Di Felice models produce the same results when the porosity is 0.4. • Koch and Hill model The model by Koch and Hill [48] was proposed on the basis of Lattice-Boltzmann simulations of fluid flows through arrays of spheres. The resulting equation for the D term is:   D ¼ 18ag F 0 þ 0:5 F 1 ag Rep

ð39Þ

rffiffiffiffiffi 8 as 135 > > as ln ðas Þ þ 16:14as þ 1 þ 3 > > 64 2 < F0 ¼ 1 þ 0:681as −8:48a2s þ 8:16a3s > > > 10a s > : a3s F 1 ¼ 0:0673 þ 0:212 as þ 0:0232 a−5 g

for asolid ≤0:4

a

Diameter for the cylindrical geometry.

!

1 þ 3 ag as þ 8:4Rep −0:343 ag



1 þ 103as Rep −0:5−2as

ð43Þ • Tenneti et al. model Tenneti, Garg and Subramaniam [86] employed simulations of mono-dispersed particles through Direct Numerical Simulations (DNS). The authors stated it is more accurate in the particle Reynolds number range of 100 to 300. Its formulation yields: D ¼ 18ag F T

ð44Þ

where  0:687

 pffiffiffiffiffi 1 þ 0:15 ag Rep 5:81as 0:48 3 as 0:61a3s 3 þ þ þ a a Re 0:95 þ g p s a4g a3g a3g a2s

ð40Þ 3. Simulation procedure

for asolid N0:4 ð41Þ

We performed the simulations with the commercial CFD program Ansys Fluent 19.2. We developed the computational geometries of the

Table 6 Numerical methods and under relaxation factors (URF) employed in the simulations.

Table 4 Dimensions of the simulated geometries.

Inlet widtha Bottom widtha Column widtha Column depth Base angle Height of the inclined part Total height

 pffiffiffiffiffi as þ a2g 1 þ 1:5 as þ a2g

0:413 Rep 24 a2g

ð45Þ

• Beetstra et al. model

Dimension

F B ¼ 10

FT ¼

where

ð42Þ

Value Pseudo-2D

Cylindrical

0.9 cm 1.5 cm 15.2 cm 1.5 cm 60° 11.9 cm 60 cm

3 cm 6.3 cm 15.2 cm – 45° 10.8 cm 55.8 cm

Solver

Pressure-based

Scheme Gradient Pressure Momentum Volume fraction Turbulence Time Pressure URF Density URF Body forces URF Momentum URF Volume fraction URF Turbulent kinetic energy URF Specific dissipation rate URF Discrete phase sources URF

Phase-coupled SIMPLE Least Squares Cell Based PRESTO! QUICK QUICK First order upwind First order upwind 0.3 1 1 0.2 0.99 0.8 0.8 0.5

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

6

F. Marchelli et al. / Powder Technology xxx (2019) xxx

Fig. 3. Distribution of the particle Reynolds numbers for upwards-travelling particles in the pseudo-2D (left) and cylindrical (right) spouted beds.

two SBs exactly matching those of the experimental devices [67,69]. Fig. 2 shows the meshes of the two geometries, while Table 4 summarises their dimensions. The first mesh (pictures on the left) has 3480 cells and 4972 nodes, while the second mesh (pictures on the right) has 13,696 cells and 14,921 nodes. For the pseudo-2D SB, the ratio is about 4 in the parallelepipedal part, while it decreases along the pyramidal part due to the converging geometry. For the cylindrical SB, the ratio is 2.5 in the vertical direction, while the value in the radial direction varies in the different parts of the SB due to the features of the structured mesh. The average ratio is 2.1, calculated considering the average volume of a cell. The grids were chosen through preliminary simulations that guaranteed the mesh-independence of the results. With regard to the pseudo-2D SB, these simulations are detailed in the previous publication [57], and are hence not reported again here. Regarding the cylindrical SB, the details of the simulations are provided in the Supporting information. The simulations begun with the particles settled in the lower part of the SB, with different initial heights. Air starts flowing from the inlet at a fixed velocity, passes through the particles exchanging momentum with them, and leaves the unit from its top, where the pressure is set as atmospheric. Table 5 summarises the properties and parameters of the simulated particles. These can have a significant influence on the overall fluid dynamics of SBs, as we showed in our previous sensitivity analysis of the spouting of glass particles in a pseudo-2D device [57]. Since all the present simulations also reproduced glass particles, we retained the same values for the particle-related parameters. Preliminarily tuning them to achieve a perfect match of the experimental data would have been pointless for the scope of this work, as the best performing set of parameters is most likely also depending on the employed drag model. We also decided not to set different values for particle-particle and particle-wall interactions, which was the most neutral choice: since we are reproducing data from the literature, it was not possible to know how exactly the parameters may differ. Besides, the involved wall materials (iron, glass and poly (methyl methacrylate)) are not expected to produce collisional phenomena that markedly differ from those happening between particles. The CFD and DEM time steps were chosen in order to satisfy common heuristics [87]. For the simulations of the pseudo-2D unit, they were respectively 5 × 10−5 s and 1 × 10−5 s. For the simulations of the cylindrical unit, we increased them to 1 × 10−4 s and 5 × 10−5 s due to the larger sizes of the particles and cells. The chosen numerical methods and under-relaxation factors [88] are displayed in Table 6. To reduce the dependency of CFD-DEM simulations on the ratio between cells and particles' sizes, Fluent can employ a node-based

averaging algorithm for the discrete phase quantities. This way, the effects of a particle's presence are not limited to the cell that contains it, but also affect the nearby cells. The equation is: φnode ¼

X


 ! ! w x i − x node φp

ð46Þ

i

where φp is the particle variable and φnode is its accumulation on the ! ! node for all of the parcels. The product wð x i − x node Þ is the weighting function or kernel, calculated in this Gaussian form: ! ! 2 !
 h3=2  x i − x node  ! ! w x i − x node ¼ exp −h π Δx2

ð47Þ

In this, Δx is the length scale of the cell containing the parcel, while h is a parameter that controls the width of the distribution. In our previous work [57], we showed that simulations are quite insensitive to the value of h, but the use of this algorithm is of key importance for the simulations to converge. In this work, the default value (1) was retained for the parameter h. The comparison is based on experimental data reported in the experimental articles [66–69]. First, we assessed the accuracy of the drag models in predicting the formation of the fountain for the reported value of the minimum spouting velocity, which is a key parameter in the operation of SBs. Then, we compared the results obtained through the simulations to the available experimental data (profiles of particle velocity and porosity). We wish to emphasise that the literature does not provide many available experimental data of SBs that are also suitable to be simulated through CFD-DEM models. This is because some cases feature too many particles, while for others it is hard to achieve an acceptable ratio between the sizes of the particles and of the cells. In fact, some researchers validated their SB models reproducing the gas velocity – bed pressure drop plot, which however was too timeconsuming to perform for seven different drag models. The comparison would also have been enhanced by assessing the variation of more variables, such as the root mean square of the velocity or the pressure fluctuations. Unfortunately, these were not available, so we focused on the variables that provided a more direct confirmation of the accuracy of each model in reproducing the spouting regime. To make an adequate comparison with the experimental profiles, we performed a time-averaging of our data. First, we started the simulations with the particles in a packed state, with the selected initial bed height and inlet gas velocity. We let the simulation run for 1 s (2 s for

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

F. Marchelliet al. / Powder Technology xxx (2019)xx x

7

Fig. 4. Plots of the dimensionlessdrag function D (on the left) and its normalisedvalue (on the right) as a function of the particleRe ynolds number and porosity for each drag model.

Pleasecite this articlea s: F. Marchelli,Q . Hou, B. Bosio, et al., Comparisonof different drag modelsin CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

8

F. Marchelli et al. / Powder Technology xxx (2019) xxx

Table 7 Formation (or lack thereof) the spouting regime with the inlet gas velocity equal to ums and 1.3ums for the various studied setups and drag models. Model

Wen-Yu Gidaspow Di Felice Rong Koch-Hill Beetstra Tenneti

7.5 cm bed, pseudo-2D

10 cm bed, pseudo-2D

20 cm bed, cylindrical

ums

1.3ums

ums

1.3ums

ums

1.3ums

No Yes No No No No No

Yes Yes No No Yes Yes Yes

No Yes No No No No No

Yes Yes No No Yes Yes No

No No No No No No No

Yes Yes No No Yes Yes No

Table 8 Operative conditions of the simulations of the pseudo-2D spouted bed. Bed height (cm)

uinlet (m/s)

u/ums

10

26.68 29.89 32.93 22.63 28.88

1.74 1.94 2.14 1.84 2.08

7.5

the cylindrical SB) to let the spouting regime establish. Then, we enabled the time-averaging algorithm and let the simulation run for another second. Preliminary simulations showed that a further increase of the time-averaging period did not affect the obtained data [57]. In light of the results of the previous sensitivity analysis [57], we wish to emphasise that CFD-DEM simulations of SBs display a remarkable dependency on the numerous sub-models and parameters. Hence, the accuracy of the results is influenced by the overall set-up, and it may be misleading to declare that a drag model is the best one, as a change in other sub-models may produce better results with another drag model. 4. Results and discussion 4.1. Magnitude of the drag force Before analysing the results of the simulations, it is good to understand the role that the drag models play. As the equations showcase, the various force definitions only differ in the dependency on the

porosity and the particle Reynolds number. The solid (and hence gas) volume fraction can be observed in Fig. 1; its ranges in the spout, fountain and annulus are quite constant for all SBs. Typically, the central spout and the fountain are quite dilute, with solid volume fraction values up to 0.3, while the spout is comparable to a packed bed, with solid volume fraction values of about 0.6. The involved particle Reynolds numbers depend on the employed particles and the chosen gas velocity. Since SBs commonly operate with coarse Geldart-D particles, high gas velocities are needed for their fluidisation. The combination of high diameters and velocity gradients leads to high values of the particle Reynolds number. This is especially true in the bottom of the spout, where the great velocity gradient between the two phases creates the vertical acceleration of particles. Conversely, in the fountain the gradient is generally lower, as the gas flow gets dispersed in the radial direction. Fig. 3 shows the distribution of the particle Reynolds number for the two studied geometries, obtained from preliminary simulations, evaluated in steady spouting conditions. The plots only show Rep for the particles travelling upwards, i.e. those in the spout and in the most internal part of the fountain, which is where the drag force has the strongest effect. The distribution for the particles travelling downwards (not pictured) features significantly lower values. These plots clearly show that the operative Rep values are quite high: for the pseudo-2D SB most of the particles have Reynolds numbers around 1000, while those in the cylindrical SB have a wider Rep distribution, with the majority of them exceeding 1000. The higher magnitude is caused by the coarser particle diameter (4 mm against 2 mm) and the higher inlet gas velocity, while the wider distribution is likely due to the higher bed height. The information about the porosity and Reynolds number ranges is useful to understand where and how the differences among the drag models will play the most relevant roles. Fig. 4 shows a comparison of the dimensionless D term for some operative conditions. In each case, the plot on the left is based on the magnitude of the term D for each model, while on the right it is normalised based on the value of the Tenneti drag model. This choice is because the Tenneti drag force has, averagely, the lowest magnitude. From Fig. 4, in all cases the drag force increases with the particle Reynolds number and decreases with the porosity. The Tenneti model always produces the lowest force magnitude, while the highest magnitude is obtained with either the Gidaspow or the Koch-Hill models. Both models also have quite evident steps when they switch from one equation to another. It is also interesting to note that the Wen-Yu model (and hence also the Gidaspow model) produces low drag forces when the

Fig. 5. Configuration of the particles (coloured by their vertical velocity) after 2 s of simulation (top) and time-averaged contour of the particle volume fraction (bottom) with u/ums = 1.74 and Hb = 10 cm in the pseudo-2D spouted bed.

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

F. Marchelli et al. / Powder Technology xxx (2019) xxx

9

Fig. 6. Time-averaged axial profile of the particles' vertical velocity (Hb = 10 cm, u/ums = 1.74, pseudo-2D spouted bed).

Fig. 8. Time-averaged vertical profile of the particles' vertical velocity (Hb = 10 cm, u/ums = 1.94, pseudo-2D spouted bed).

porosity is high. This is the condition in which the model is considered most accurate [73]. The Di Felice and Rong drag models yield similar values for the drag force, identical when the porosity approaches 0.4.

expected that the different drag definitions will produce different outcomes at the experimental minimum spouting conditions. Hence, we first ran simulations setting the inlet gas velocity to ums as reported in the original papers [67,68]. Due to the hysteretic behaviour of SBs, we first let the simulation run for 1.5 s with a velocity value equal to 1.3ums. Then, we decreased the velocity to the reported value of ums. Table 7 summarises the results: ‘yes’ means that at the end of the simulation the particles were in the spouting regime, while ‘no’ means that they were in a fixed-bed state. Surprisingly, given its simplistic nature, the Gidaspow model correctly predicts the particles' fluidisation in all cases, besides in the cylindrical SB at ums. However, the authors of the original paper did not report the actual ums values, but rather just stated that they fit the Mathur-Gishler equation [1] with an average relative error of 15% [68]. This can justify the lack of fluidisation for that single case, confirming the adequacy of the Gidaspow drag force in these conditions.

4.2. Prediction of the minimum spouting velocity In SBs, the minimum spouting velocity (ums) is the lowest gas velocity value for which the fountain does not collapse. Due to the hysteretic behaviour of SBs, it is usually slightly lower than the onset velocity (i.e. the minimum value for which the fountain establishes). It is a key parameter in the operation of SBs that needs to be known before running the device, and is affected by several variables (geometric features, bed height, particle size, density and shape). In a previous work, we showed that CFD-DEM simulations can reasonably respond to the variation of these variables, but the exact prediction of the operational gas velocities is not always achieved [89]. Given their different magnitudes, it can be

Fig. 7. Time-averaged radial profile of the particles' vertical velocity in the spout (left) and annulus (right) (Hb = 10 cm, u/ums = 1.74, pseudo-2D spouted bed).

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

10

F. Marchelli et al. / Powder Technology xxx (2019) xxx

Fig. 9. Configuration of the particles (coloured by their vertical velocity) after 2 s of simulation (top) and time-averaged contour of the particle volume fraction (bottom) with u/ums = 2.14 and Hb = 10 cm in the pseudo-2D spouted bed.

The Wen-Yu, Koch-Hill and Beetstra models behave similarly: they do not create the spouting regime at the minimum spouting conditions, but they behave correctly with the slightly increased velocity. Finally, the Di Felice, Rong and Tenneti models do not predict the fluidisation in both conditions.

Even from this preliminary analysis, there seems to be a clear trend that the models behave similarly in different geometries and bed heights. The behaviours are coherent with the theoretical magnitudes of the different drag forces, which can be inspected through plots of Fig. 4 with ag equal to 0.4 (i.e. the initial packed bed conditions).

Fig. 10. Time-averaged vertical (left) and radial (right) profiles of the particles' vertical (Hb = 10 cm, u/ums = 2.14, pseudo-2D spouted bed).

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

F. Marchelli et al. / Powder Technology xxx (2019) xxx

11

Fig. 11. Time-averaged vertical profiles of the particles' vertical velocity (Hb = 7.5 cm, u/ums = 1.84 (left), 2.08 (right), pseudo-2D spouted bed).

When Rep exceeds 1000, the Gidaspow model yields the highest force magnitude, followed by Wen-Yu, Koch-Hill and Beetstra. When Rep is equal to 1000, Gidaspow continues to yield the highest force magnitude up to a porosity value of 0.8, i.e. when the equation shifts to the Wen-Yu formulation. The Wen-Yu, Koch-Hill and Beetstra models yield drag functions of the same order of magnitude, but their slightly lower values with respect to the Gidaspow model are enough to impede the initiation of the spouting regime. This confirms the great reliability of the Ergun equation. The Rong, Di Felice and Tenneti models yield significantly lower values. Rong and Di Felice produce a lower force than the Tenneti model only in a small range at low Rep values. This is likely why Tenneti model correctly produces the spouting regime for the 7.5 cm bed when u = 1.3ums. Indeed, given the lower initial bed height, a lower gas velocity is necessary, and hence Rep is averagely lower. Based on these results, we preliminarily conclude that the Gidaspow model is the most suitable when simulating a SB at minimum spouting conditions. All the other models yield underpredicted values for the drag force in these conditions. 4.3. Pseudo-2D spouted bed For the pseudo-2D SB, we compared the results with the experimental vertical and radial profiles of the vertical velocity of particles (Vz). The vertical profile is referred to the central axis of the SB, while the radial profile is calculated at a height of 9.12 cm from the bottom of the bed. Table 8 summarises the operative conditions of the simulations. 4.3.1. Initial bed height: 10 cm This Section presents the results of the simulations performed with a 10 cm initial bed height, in the operative conditions described in Table 8. For the lowest and highest gas velocities, a visual depiction of the results is also presented. For the inlet gas velocity of 26.68 m/s, Fig. 5 presents instantaneous snapshots of the particle configuration after 2 s of simulation and timeaveraged contours of the particle volume fraction. From a qualitative point of view, all the models seem valid: the three typical zones of SBs (annulus, spout and fountain) can clearly be distinguished in all cases. By comparing the configurations of the particles with the experimental

snapshot, the Di Felice, Rong and Tenneti models result in a more comparable spouting regime, with a very low fountain and a denser spout. The other models produce significantly higher fountains. In all cases, at the top of the fountain an accumulation zone forms, with an associated porosity of about 0.8. To analyse the results in better detail, Fig. 6 provides the timeaveraged axial profiles of the particles' vertical velocity, compared with the experimental values. From the figure, it is apparent that none of the models can give a perfect prediction of the particle velocity profiles. Nonetheless, the Rong model seems to represent the best compromise considering the maximum velocity, profile shape and fountain height. The errors are −9.8% for the peak velocity, and 14.8% for the fountain height. The Di Felice model creates a similar fountain height, but averagely lower velocity values (−15% error). This is consistent with the models' obtained D values (as per Fig. 4): they are equivalent when the porosity is 0.4 [59], but the Rong model predicts slightly higher drag force values for higher porosity values. In a SB, the drag force is prominent in the spout, where the porosity is high, which justifies the observed profiles. The Tenneti model produces a slightly lower maximum Vz, but a higher fountain than those of the Di Felice and Rong models. This is coherent with the fact that, for low Reynolds number, the Tenneti model results in higher drag force values than the Di Felice and Rong models. This is observable in the fourth plot of Fig. 4; the difference is minimal, but enough to cause a difference. The Tenneti model also creates the worst profile shape in the first part of the spout, damping the abrupt velocity increase of this region. Again, this is coherent considering that the Tenneti model produces the weakest force in most cases. The Wen-Yu model gives a slightly underpredicted maximum Vz, and a largely overpredicted fountain height, although both are lower than those of the Gidaspow model. The introduction of the Ergun equation that the Gidaspow model operates is apparently counterproductive in these conditions. This confirms that, although the drag is prominent in the spout and fountain, its effect in the more packed zones also has a marked influence on the results. In fact, if this effect were negligible, the Wen-Yu and Gidaspow models would produce identical results. Finally, both the Koch-Hill and Beetstra models provide largely overpredicted profiles. The plots of Fig. 4 justify these trends, as these two models produce the most intense force magnitudes when the

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

12

F. Marchelli et al. / Powder Technology xxx (2019) xxx

Fig. 12. Central section of the configuration of the particles (coloured by their volume fraction) after 3 s of simulation (top) and time-averaged contour of the particle volume fraction (bottom) with u/ums = 1.3 and Hb = 20 cm in the cylindrical spouted bed.

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

F. Marchelli et al. / Powder Technology xxx (2019) xxx

13

Fig. 13. Axial (left) and radial (right) profiles of the time-averaged solid volume fraction in the cylindrical spouted bed (u/ums = 1.3).

porosity and Reynolds number are both high. It is also worth noting that the Di Felice and Rong force magnitudes are higher than those of WenYu (and hence Gidaspow) in high-porosity conditions, as clearly observable in the last plot of Fig. 4. Nevertheless, Di Felice and Rong result in lower fountain heights. This may indicate that this parameter depends more on the acceleration that the particles undergo in the spout. Indeed, the gas travels more slowly when it reaches the fountain area, not generating intense drag effects anymore. Fig. 7 shows that the results of the radial velocity profiles feature a lower variability. Overall, the Di Felice and Rong models are confirmed as the best performing. Most notably, all the models predict very well the spout width, which is commonly considered as the point where the vertical velocity of the particles switches from negative to positive. This, interestingly, proves that the various definitions of D have little effect on the spout width. This partially contrasts with the previous observation that the drag force also has an influence on the packed zone. We conclude that this influence is likely more on the velocity of the particles and the spout-annulus mass exchange rate. In addition to a direct drag effect, the more particles are dragged upwards in the spout, the faster the particles in the annulus will move downwards to replace them. Hence, one would impulsively assume that the models that yield the most intense D values in the spout (Koch-Hill and Beestra) and the highest upwards velocities, should also yield the lowest negative velocities in the annulus. This is however not the case (as per the left plot of Fig. 7). The causes may be twofold. First, it may mean that the drag effect not only moves particles upwards in the spout, but also slows down their descending movement in the annulus. Second, the velocity of particles in the annulus must not only be affected by the mean velocity of particles in the spout, but also by the spout's porosity and width. For example, the Di Felice model yields a moderate particle velocity in the annulus, but a slightly larger spout than the Rong model. As a result, in the simulations with the Di Felice model the particles descend slightly more quickly in the annulus. Upon increasing the gas velocity to 29.89 m/s, the relative behaviour of the models changes slightly: Fig. 8 showcases the vertical profile of Vz under this condition. The horizontal experimental profiles were not available for this case. As in the previous case, the Di Felice model gives the best prediction of the fountain height (10% error), although

it still underestimates the maximum velocity (−18% error). Conversely, the Rong model's accuracy worsens, as it starts displaying a saddleshape that is not encountered in the experimental data. Tenneti and Wen-Yu are the only other models that predict a flat-enough shape. However, the former yields again a very low velocity in the first part of the spout, while the latter largely overestimates the fountain height. The other models give even worse results, with very evident saddle shapes and even more overestimated fountain heights. The obtained relative trends are similar to those of the previous case, and again coherent with the magnitude of the D functions. The Tenneti profile lies again lower than the Rong and Di Felice profiles in the first part of the spout, and then ends very close to them, more than in the previous case. This difference is due to the averagely higher Reynolds number. Moreover, the Beetstra profile gets higher than that of the Gidaspow one, which in turn becomes more similar to that of Wen-Yu. This can again be explained considering that, in these simulations, the yielded D values are shifted more to the left along the Rep-D plots of Fig. 4. Finally, we increased the inlet gas velocity to 32.93 m/s, reaching a u/ ums ratio of 2.14. Fig. 9 provides the visual depiction of the spouting regime and the time-averaged volume fraction. Upon a qualitative comparison with the experimental particles' configuration, the Rong and Tenneti models produce the most realistic particle configurations. With the Di Felice model, the fountain is smaller in the horizontal direction. The four remaining models largely overpredict the bed expansion. For a more quantitative comparison, Fig. 10 shows the timeaveraged velocity profiles. The vertical profile permits to draw very similar conclusions to the previous case: The Di Felice model underestimates the maximum velocity (−15.3% error), but gives the best prediction of the fountain height (9.8% error) and seems the overall best compromise. The Rong and Tenneti models create good results in the first and second part of the spout respectively, while they behave similarly in the fountain, overestimating its maximum height by about 16%. All the other models continue to be inadequate for these highvelocity working conditions. The effect of the averagely-higher Rep observed in the previous case affect these results in a similar way: the average velocity obtained through the Beetstra model is even higher than that of the Gidaspow model, the Gidaspow and Wen-Yu model yield nearly identical profiles in the second part, and the Tenneti profile

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

14

F. Marchelli et al. / Powder Technology xxx (2019) xxx

Fig. 14. Axial (left) and radial (right) profiles of the time-averaged particles' vertical velocity in the cylindrical spouted bed (u/ums = 1.3).

does not end beyond the Rong model anymore. All these effects are very consistent with the plots of Fig. 4. A more cryptic consequence is the formation of a saddle-shaped profile, which happens with all the models besides Di Felice, Tenneti and Wen-Yu. Such a decrease in the velocity in proximity of the end of the spout may be caused by an accumulation of particles. These three models are characterised by very different drag functions, so it is quite difficult to understand what makes them behave similarly. It is also quite intriguing that the Rong model produces the saddle shape while the Di Felice model does not, despite their very similar formulations. As in the first case, the radial profile is less affected by the choice of the drag model, but some differences can be observed. Tenneti gives the best results in the spout, but slightly overestimates its width. None of the model can produce a good matching of the experimental profile in the annulus, but all the results lie quite close to it. 4.3.2. Initial bed height: 7.5 cm We ran further simulations with an initial bed height of 7.5 cm and u/ums ratios of 1.84 and 2.08 to assess whether similar trends could be observed even with a lower bed, and hence averagely lower Rep. We chose these two u/ums values because, among the reported experimental data, they were the closest to the ones employed for the 10 cm bed (1.74 and 2.14). Fig. 11 summarises the results. Indeed, the different drag models produce results akin to the previous cases. Some differences can be observed: the discrepancy between the Di Felice and Rong model is less evident, which is coherent with the fact that the average particle Reynolds number is lower. Again, the Rong and Di Felice models appear more suitable for the moderate u/ ums ratio: the errors on the peak velocity and fountain height are respectively −10.5% and 23% with the Rong model and −17.3% and 15.5% with the Di Felice model. Interestingly, both models generate a peak in the velocity at the beginning of the spout, while the other models generate it in the fountain.

With the high u/ums ratio, the Di Felice model respectively yields errors on the peak velocity and fountain height of −5.5% and 26.6%, whereas those of the Rong model are 2.3% and 39.3%. In this case, the Tenneti and Wen-Yu models also produce very similar and still reasonable results: their errors on the peak velocity are similar to that obtained with the Di Felice model, but the error on the fountain height is higher. Finally, as in the previous cases, the remaining models give largely overpredicted velocity profiles, confirming their inadequacy for the present working conditions. For this u/ums ratio, all the models create an unrealistic saddle shape in the profile, which was not so evident in the previous cases. However, the behaviour of the Di Felice and Rong models observed at the low u/ums value is not encountered anymore: all models feature the peak velocity in proximity of the fountain. As the simulations at u = ums seemed to suggest, these results confirm that the behaviour of the drag models is not significantly affected by the bed height. The choice should be primarily based on the u/ums ratio. The generally lower accuracy of the results for the lower bed may also be because, in this case, the majority of particles lie in the zone where their size approaches that of the mesh cells. 4.4. Cylindrical spouted bed This Section presents the results of the simulations of the cylindrical SB by Olazar and colleagues [68,69]. As we wrote in Section 4.2, some models did not generate the spouting regime in the reported experimental conditions, so the comparison is only based on the models that did: Wen-Yu, Gidaspow, KochHill, Beetstra. The other three models generate lower drag forces in packed conditions and therefore they could not establish the spouting regime for u = 1.3u ms . We performed the simulations with an inlet gas velocity of 45.40 m/s, i.e. with a u/u ms ratio equal to 1.3. Given the lack of experimental data at 1.3u ms , the comparison is in some cases performed with the data reported at

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

F. Marchelli et al. / Powder Technology xxx (2019) xxx

1.02ums, for which the trend and values are quite similar. In this case, the available experimental data regarded the particles' vertical velocity (Vz) and volume fraction (as). The vertical profiles are referred to the central axis of the SB, while the radial profiles are calculated at a height of 8 cm from the bottom of the bed. First, Fig. 12 gives a visual representation of the spouting regime through the four suitable models; above is the instantaneous depiction of the particles' configuration at the end of the simulations, while below is the time-averaged volume fraction contour. From a qualitative point of view, the four models create a reasonable fluid dynamic regime, and the distinguished regions of the SB are clearly visible. A more accurate inspection of the time-averaged volume fraction contours highlights that the Gidaspow and Koch-Hill models generate higher and similar fountains. Conversely, the fountains of the Wen-Yu and Beetstra models are lower, and with the latter model also less porous. Considering the obtained magnitude of the drag functions, this seems to reinforce that, despite the high Reynolds numbers of the ascending particles, the overall behaviour is greatly affected by the packed zone. If this were not the case, the Beetstra model would produce a higher fountain than the Gidaspow and Wen-Yu models, as its D is higher at high Rep and ag. By contrast, its yielded D is lower than those of the Wen-Yu and Gidaspow models in packed conditions. The Koch-Hill model results in high D values both in packed conditions and at moderate Rep (as in the annulus) and in dilute conditions and at large Rep (as in the spout and fountain). Consequently, it makes particles travel higher and, as the next plots will show, faster. Fig. 13 provides a visual depiction of the solid volume fraction profiles. This variable is least affected in the radial direction: the models produce nearly identical profiles, all very close to the experimental data. The prediction of the vertical profile is less accurate, although all the models create qualitatively reasonable trends. In this case, the Beetstra model is clearly the best performing one, as it yields the closest profile to the experimental data, especially for the peak at the top of the spout and for the fountain height. The other three models further underestimate the former and overestimate the latter. The trends are similar for all models, with an absolute peak at the upper end of the spout, and another relative peak in the centre of the fountain. The velocity profiles, depicted in Fig. 14, permit a further assessment of the models' performances. As it is apparent from the vertical profile, all of them generate a peak in the velocity that is axially higher than in the experiments, and with a significantly lower peak value. However, the Koch-Hill and Beetstra model can at least generate a plausible profile, with a very pronounced peak at the beginning of the spout, followed by a swift decline of the velocity. Conversely, the Wen-Yu and Gidaspow models produce nearly flat profiles, like those of the pseudo-2D SB. This is an unwanted result, which should discourage their use in the present conditions. The previous reasonings are coherent with the behaviours of the Beetstra and Koch-Hill models. Clearly, Beetstra produces such an intense velocity peak near the inlet thanks to its high magnitude at high Rep and ag. The decrease in velocity and accumulation of particles at the top end of the spout foster each other, resulting in the lowest fountain. Assessing the reason for the flatter profiles of the Gidaspow and Wen-Yu models is harder, and most likely cannot be ascribed to a single cause. Again, the marked difference between Gidaspow and Wen-Yu confirms the considerable role of the drag force in the packed zone of the SB. The radial profiles show fewer variations, especially in the annulus (bottom plot). The models behave quite similarly: they all predict a 1cm-larger spout and underestimate the velocity magnitude in both the upwards and downwards direction. This is reasonable, as the faster particles move upwards in the spout, the faster they will also move downwards in the annulus, to replenish the bottom of the bed. The fact that the simulations overpredict the spout width probably also causes the

15

underestimation of the particles' velocity in the spout, as the gas flow is distributed over a larger section. Considering these findings, we conclude that no model is entirely accurate for the simulation of this SB. Indeed, all of them underestimated the drag magnitude, especially in the lowest part of the spout. Nonetheless, the Beetstra model yielded reasonable predictions for both the solid volume fraction and vertical velocity and hence can be considered the best compromise.

5. Conclusions In this work, we employed CFD-DEM simulations to study the behaviour of seven different drag models: three classic models (Wen-Yu, Gidaspow, Di Felice), three developed with Lattice-Boltzmann simulations (Rong, Beetstra, Koch-Hill), and one developed with Direct Numerical Simulations (Tenneti). We simulated two different spouted beds whose data are reported in the literature, one pseudo-2D and one cylindrical, containing Geldart-D particles. For the studied cases, the simulations highlighted that neither the column depth nor the bed height affects the adequacy of these models. Conversely, the u/ums ratio seems to be the key variable, which is reasonable considering its link to the particle Reynolds number, one of the key variables of the drag force. In all the studied cases, the Gidaspow model was the only one able to predict the particles' fluidisation for the reported ums. For a 30% increased gas velocity value, the Wen-Yu, Beetstra and Koch-Hill model correctly produced the spouting regime too. In the cylindrical spouted bed, for which the available experimental data were in the range 1–1.3ums, the Beetstra model was deemed as the best compromise to reproduce the volume fraction and velocity profiles. However, all the studied models underestimated the peak velocity and overestimated the fountain height. For the pseudo-2D spouted bed, for which the available data were in the range 1.74–2.14ums, the situation was opposite. The Wen-Yu, Gidaspow, Beetstra and Koch-Hill models largely overestimated the fountain height and, in some cases, the peak velocity. The Di Felice, Rong and Tenneti models were instead able to generate reasonable profiles. The Rong model performs best when the u/ums ratio is 1.74, while the Di Felice model produces the best results when this value is further increased, especially above 2. Despite the very coherent results for the prediction of ums in both geometries, it cannot be ruled out that there is also an influence of wall effects in the pseudo-2D spouted bed. Nonetheless, such an effect should not significantly affect the drag force, as the low mesh resolution probably does not allow accounting for the wall effects on the gas flow. The obtained results are coherent with the theoretical force magnitudes obtained with the different models. Interestingly, there is a relevant influence of the drag force also in the packed zone, where its magnitude is at its lowest. This in turn affects the overall behaviour of the spouted beds. We wish to point out that even the models that we deemed as the best performing for each case were not able to provide a perfect reproduction of the experimental data. This may be partially fixed by modifying other submodels but may ultimately mean that the typical operative conditions of spouted beds (high Re p and porosity) are beyond the range of application of common drag models. In fact, these are typically developed employing data from fixed beds, at intermediate Rep values. We hope this work will foster more research efforts in this direction. Given the previous considerations, new models (or adequate combinations of new and old models for the different regions) should work properly both at low and high values of Re p and porosity. Nonetheless, and despite the scarcity of suitable experimental data in the literature, these findings represent valid recommendations for researchers who aim to simulate spouted beds.

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

16

F. Marchelli et al. / Powder Technology xxx (2019) xxx

Nomenclature Ap projected particle surface area [m2] aq volume fraction of phase q [−] CRL rotational lift coefficient [−] Cω rotational drag coefficient [−] D drag function [−] dij deformation tensor [−] dp particle diameter [m] eij unit vector for parcels i and j [−] F0, F1 functions in the Koch-Hill drag model [−] FB function in the Beetstra drag model [−] Fc contact force [N] Fd drag force [N] Ffriction rolling friction force [N] Fml Magnus lift force [N] Fnormal normal force [N] Fpg pressure gradient force [N] Frolling friction force [N] FT function in the Tenneti drag model [−] Fvm virtual mass force [N] floss loss factor [−] Gk production of turbulent kinetic energy [J/(m3·s)] g gravitational acceleration [m/s2] h width of the Gaussian distribution [−] I identity tensor [−] Ip particle moment of inertia [kg·m2] K spring-dashpot constant [N/m] Kpq gas-solid exchange coefficient [kg/(m3·s)] k turbulent kinetic energy [m2/s2] m mass of a particle [kg] mij reduced mass for parcels i and j [kg] np number of particles per cell [1/m3] p pressure [Pa] Rpq gas-solid momentum exchange [N/m3] Re Reynolds number [−] Rep particle Reynolds number [−] Reω rotational Reynolds number [−] t time [s] tcoll collision time scale [s] tq stress-strain tensor [Pa] Uq phase-weighted velocity [m/s] u fluid phase velocity [m/s] uij relative velocity for parcels i and j [m/s] up parcel velocity [m/s] Vz time-averaged vertical velocity of the particles [m/s] w weighting function [−] xnode position vector of a generic node [m] z vertical coordinate [m] Greek symbols β function of ag and Rep in the Di Felice and Rong models [−] γ damping coefficient [−] Δt CFD time step [s] ΔtDEM DEM time step [s] Δx length scale of a cell [m] δ parcels overlap [m] ε turbulence dissipation rate [J/(kg·s)] η restitution coefficient for the dashpot term [−] μ dynamic viscosity [Pa·s] μf friction coefficient [−] μrolling rolling friction coefficient [−] ν kinematic viscosity [m2/s] ρ density [kg/m3] φp generic particle variable φnode accumulation of a generic particle variable in a node Ω relative particle-fluid angular velocity [rad/s]

ωp

parcel rotational velocity [rad/s]

Subscripts g relative to the gas phase ms minimum spouting p particle q generic phase s relative to the solid phase t turbulent Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The authors are grateful to the Australian Research Council (IH140100035, DE180100266) for the financial support. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.powtec.2019.10.058. References [1] K.B. Mathur, P.E. Gishler, A technique for contacting gases with coarse solid particles, AICHE J. 1 (1955) 157–164, https://doi.org/10.1002/aic.690010205. [2] R.C. de Brito, T.F. de Pádua, J.T. Freire, R. Béttega, Effect of mechanical energy on the energy efficiency of spouted beds applied on drying of sorghum [Sorghum bicolor (L) moench], Chem. Eng. Process. Process Intensif. 117 (2017) 95–105, https://doi. org/10.1016/J.CEP.2017.03.021. [3] G. Nakamura Alves Vieira, F. Bentes Freire, J.T. Freire, Control of the moisture content of milk powder produced in a spouted bed dryer using a grey-box inferential controller, Dry. Technol. 33 (2015) 1920–1928, https://doi.org/10.1080/07373937. 2015.1075999. [4] J.H. van Laar, H. Bissett, J.C. Barry, I.J. van der Walt, P.L. Crouse, Deposition of SiC/Si coatings in a microwave plasma-assisted spouted bed reactor, J. Eur. Ceram. Soc. 38 (2018) 1197–1209, https://doi.org/10.1016/J.JEURCERAMSOC.2017.10.030. [5] X. Liu, W. Zhong, A. Yu, B. Xu, J. Lu, Mixing behaviors in an industrial-scale spoutfluid mixer by 3D CFD-TFM, Powder Technol. (2016) https://doi.org/10.1016/j. powtec.2016.10.046. [6] C. Beltramo, G. Rovero, G. Cavaglià, Hydrodynamic and thermal experimentation on square-based spouted beds for polymer upgrading and unit scale-up, Can. J. Chem. Eng. 87 (2009) 394–402, https://doi.org/10.1002/cjce.20172. [7] J.F. Saldarriaga, R. Aguado, A. Atxutegi, J. Bilbao, M. Olazar, Kinetic modelling of pine sawdust combustion in a conical spouted bed reactor, Fuel. 227 (2018) 256–266, https://doi.org/10.1016/j.fuel.2018.04.060. [8] J. Alvarez, B. Hooshdaran, M. Cortazar, M. Amutio, G. Lopez, F.B. Freire, M. Haghshenasfard, S.H. Hosseini, M. Olazar, Valorization of citrus wastes by fast pyrolysis in a conical spouted bed reactor, Fuel. 224 (2018) 111–120, https://doi.org/10. 1016/j.fuel.2018.03.028. [9] G. Rovero, A.P. Watkinson, A two-stage spouted bed process for autothermal pyrolysis or retorting, Fuel Process. Technol. 26 (1990) 221–238, https://doi.org/10.1016/ 0378-3820(90)90007-F. [10] D. Bove, C. Moliner, M. Curti, M. Baratieri, B. Bosio, G. Rovero, E. Arato, Preliminary tests for the thermo-chemical conversion of biomass in a spouted bed pilot plant, Can. J. Chem. Eng. (2018) https://doi.org/10.1002/cjce.23223. [11] I.K. Alghurabie, B.O. Hasan, B. Jackson, A. Kosminski, P.J. Ashman, Fluidized bed gasification of Kingston coal and marine microalgae in a spouted bed reactor, Chem. Eng. Res. Des. 91 (2013) 1614–1624, https://doi.org/10.1016/j.cherd.2013.04.024. [12] A. Niksiar, B. Nasernejad, Activated carbon preparation from pistachio shell pyrolysis and gasification in a spouted bed reactor, Biomass Bioenergy 106 (2017) 43–50, https://doi.org/10.1016/j.biombioe.2017.08.017. [13] Z. Wang, C.J. Lim, J.R. Grace, Biomass torrefaction in a slot-rectangular spouted bed reactor, Particuology. 42 (2019) 154–162, https://doi.org/10.1016/J.PARTIC.2018. 02.002. [14] D. Liu, E.P.L. Roberts, A.D. Martin, S.M. Holmes, N.W. Brown, A.K. Campen, N. de las Heras, Electrochemical regeneration of a graphite adsorbent loaded with acid violet 17 in a spouted bed reactor, Chem. Eng. J. 304 (2016) 1–9, https://doi.org/10.1016/j. cej.2016.06.070. [15] S. Yang, Y. Sun, J. Wang, A. Cahyadi, J.W. Chew, Influence of operating parameters and flow regime on solid dispersion behavior in a gas-solid spout-fluid bed, Chem. Eng. Sci. 142 (2016) 112–125, https://doi.org/10.1016/j.ces.2015.11.038.

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

F. Marchelli et al. / Powder Technology xxx (2019) xxx [16] L. Spreutels, J. Chaouki, F. Bertrand, B. Haut, R. Legros, Gas residence time distribution in a conical spouted bed, Powder Technol. 290 (2016) 62–71, https://doi.org/ 10.1016/j.powtec.2015.07.033. [17] A. Pablos, R. Aguado, M. Tellabide, H. Altzibar, F.B. Freire, J. Bilbao, M. Olazar, A new fountain confinement device for fluidizing fine and ultrafine sands in conical spouted beds, Powder Technol. 328 (2018) 38–46, https://doi.org/10.1016/j. powtec.2017.12.090. [18] P. Kieckhefen, S. Pietsch, M. Höfert, M. Schönherr, S. Heinrich, F. Kleine Jäger, Influence of gas inflow modelling on CFD-DEM simulations of three-dimensional prismatic spouted beds, Powder Technol. 329 (2018) 167–180, https://doi.org/10. 1016/J.POWTEC.2018.01.048. [19] T. Madhiyanon, S. Soponronnarit, W. Tia, Industrial-scale prototype of continuous spouted bed paddy Dryer, Dry. Technol. 19 (2001) 207–216, https://doi.org/10. 1081/DRT-100001362. [20] M. Jacob, ProCell technology: Modelling and application, Powder Technol. 189 (2009) 332–342, https://doi.org/10.1016/j.powtec.2008.04.035. [21] X. Bao, W. Du, J. Xu, An overview on the recent advances in computational fluid dynamics simulation of spouted beds, Can. J. Chem. Eng. 91 (2013) 1822–1836, https://doi.org/10.1002/cjce.21917. [22] X. Bao, W. Du, J. Xu, Computational fluid dynamic modeling of spouted beds, in: N. Epstein, J.R. Grace (Eds.), Spouted Spout-Fluid Beds, Cambridge University Press, Cambridge 2010, pp. 57–81, https://doi.org/10.1017/CBO9780511777936.005. [23] C. Moliner, F. Marchelli, B. Bosio, E. Arato, Modelling of spouted and spout-fluid beds: key for their successful scale up, Energies. 10 (2017) 38, https://doi.org/10. 3390/en10111729. [24] Z. Rahimi-Ahar, M.S. Hatamipour, Hydrodynamics, numerical study and application of spouted bed, Rev. Chem. Eng. (2018) 1–24, https://doi.org/10.1515/revce-20170036. [25] S. Ergun, Fluid flow through packed columns, Chem. Eng. Prog. 48 (1952) 89–94. [26] C.Y.Y. Wen, Y.H.H. Yu, Mechanics of fluidization, Chem. Eng. Prog. Symp. Ser. 162 (1966) 100–111. [27] R. Di Felice, The voidage function for fluid-particle interaction systems, Int. J. Multiphase Flow 20 (1994) 153–159, https://doi.org/10.1016/0301-9322(94)90011-6. [28] L. Li, J. Remmelgas, B.G.M. van Wachem, C. von Corswant, S. Folestad, M. Johansson, A. Rasmuson, Effect of drag models on residence time distributions of particles in a Wurster fluidized bed: a DEM-CFD study, KONA Powder Part J. 33 (2016) 264–277, https://doi.org/10.14356/kona.2016008. [29] L. Lu, A. Konan, S. Benyahia, Influence of grid resolution, parcel size and drag models on bubbling fluidized bed simulation, Chem. Eng. J. 326 (2017) 627–639, https://doi. org/10.1016/J.CEJ.2017.06.002. [30] R. Stanly, G. Shoev, Detailed analysis of recent drag models using multiple cases of mono-disperse fluidized beds with Geldart-B and Geldart-D particles, Chem. Eng. Sci. 188 (2018) 132–149, https://doi.org/10.1016/J.CES.2018.05.030. [31] Y. Wu, D. Liu, J. Ma, X. Chen, Effects of gas-solid drag model on Eulerian-Eulerian CFD simulation of coal combustion in a circulating fluidized bed, Powder Technol. 324 (2018) 48–61, https://doi.org/10.1016/j.powtec.2017.10.013. [32] N.V. Koralkar, M. Bose, Performance of drag models for simulation of fluidized beds with Geldart D particles, Adv. Powder Technol. 27 (2016) 2377–2398, https://doi. org/10.1016/j.apt.2016.11.008. [33] N. Mahinpey, F. Vejahati, N. Ellis, CFD simulation of gas—Solid bubbling fluidized bed: An extensive assessment of drag models, Comput. Methods Multiph. Flow IV, WIT Press, Southampton, UK 2007, pp. 51–60, https://doi.org/10.2495/MPF070061. [34] C. Loha, H. Chattopadhyay, P.K. Chatterjee, Assessment of drag models in simulating bubbling fluidized bed hydrodynamics, Chem. Eng. Sci. 75 (2012) 400–407, https:// doi.org/10.1016/j.ces.2012.03.044. [35] J. Xie, W. Zhong, A. Yu, MP-PIC modeling of CFB risers with homogeneous and heterogeneous drag models, Adv. Powder Technol. 29 (2018) 2859–2871, https://doi. org/10.1016/J.APT.2018.08.007. [36] J.F. Richardson, W.N. Zaki, Sedimentation and fluidisation. Part 1, Trans. Inst. Chem. Eng. 32 (1954) 35–53. [37] D. Gidaspow, Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions, Academic Press, 1994. [38] M. Syamlal, T. O'Brien, Computer simulation of bubbles in a fluidized bed, AIChE Symp. Ser. 85 (1989) 22–31. [39] H. Arastoopour, P. Pakdel, M. Adewumi, Hydrodynamic analysis of dilute gas—solids flow in a vertical pipe, Powder Technol. 62 (1990) 163–170, https://doi.org/10. 1016/0032-5910(90)80080-I. [40] W. Du, X. Bao, J. Xu, W. Wei, Computational fluid dynamics (CFD) modeling of spouted bed: assessment of drag coefficient correlations, Chem. Eng. Sci. 61 (2006) 1401–1420, https://doi.org/10.1016/j.ces.2005.08.013. [41] J. Gao, X. Lan, Y. Fan, J. Chang, G. Wang, C. Lu, C. Xu, Hydrodynamics of gas–solid fluidized bed of disparately sized binary particles, Chem. Eng. Sci. 64 (2009) 4302–4316, https://doi.org/10.1016/j.ces.2009.07.003. [42] Y. Zhang, J.M. Reese, The drag force in two-fluid models of gas–solid flows, Chem. Eng. Sci. 58 (2003) 1641–1644, https://doi.org/10.1016/S0009-2509(02)00659-0. [43] D.A. Santos, I.J. Petri, C.R. Duarte, M.A.S. Barrozo, Study of hybrid drag models for predicting hydrodynamic behaviour in a spouted bed, Can. J. Chem. Eng. 91 (2013) 1776–1785, https://doi.org/10.1002/cjce.21866. [44] J. Prieur Du Plessis, Analytical quantification of coefficients in the Ergun equation for fluid friction in a packed bed, Transp. Porous Media 16 (1994) 189–207, https://doi. org/10.1007/BF00617551. [45] S.H. Hosseini, G. Ahmadi, M. Olazar, CFD simulation of cylindrical spouted beds by the kinetic theory of granular flow, Powder Technol. 246 (2013) 303–316, https:// doi.org/10.1016/j.powtec.2013.05.016.

17

[46] 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. [47] S. Wang, B. Shao, R. Liu, J. Zhao, Y. Liu, Y. Liu, S. Yang, Comparison of numerical simulations and experiments in conical gas-solid spouted bed, Chin. J. Chem. Eng. 23 (2015) 1579–1586, https://doi.org/10.1016/j.cjche.2015.07.006. [48] D.L. Koch, R.G. Hill, Inertial effects in suspensions and porous-media flows, Annu. Rev. Fluid Mech. 33 (2001) 619–647, https://doi.org/10.1146/annurev.fluid. 33.1.619. [49] R. Beetstra, M.A. Van Der Hoef, J.A.M. Kuipers, Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres, AICHE J. 53 (2007) 489–501, https://doi.org/10.1002/aic.11065. [50] S. Pietsch, S. Heinrich, K. Karpinski, M. Müller, M. Schönherr, F. Kleine Jäger, CFDDEM modeling of a three-dimensional prismatic spouted bed, Powder Technol. (2016) https://doi.org/10.1016/j.powtec.2016.12.046. [51] J.P. Du Plessis, J.H. Masliyah, Flow through isotropic granular porous media, Transp. Porous Media 6 (1991) 207–221, https://doi.org/10.1007/BF00208950. [52] M.A. Van Der Hoef, R. Beetstra, J.A.M. Kuipers, Lattice-Boltzmann simulations of lowReynolds-number flow past mono- and bidisperse arrays of spheres: results for the permeability and drag force, J. Fluid Mech. 528 (2005) 233–254, https://doi.org/10. 1017/S0022112004003295. [53] L. Li, B. Li, Z. Liu, Modeling of spout-fluidized beds and investigation of drag closures using OpenFOAM, Powder Technol. 305 (2017) 364–376, https://doi.org/10.1016/j. powtec.2016.10.005. [54] L. Zhou, L. Zhang, L. Bai, W. Shi, W. Li, C. Wang, R. Agarwal, Experimental study and transient CFD/DEM simulation in a fluidized bed based on different drag models, RSC Adv. 7 (2017) 12764–12774, https://doi.org/10.1039/C6RA28615A. [55] S.R. Dahl, C.M. Hrenya, Size segregation in gas–solid fluidized beds with continuous size distributions, Chem. Eng. Sci. 60 (2005) 6658–6673, https://doi.org/10.1016/j. ces.2005.05.057. [56] H. Zhang, S. Li, Study on drag force coefficients in modeling granular flows in a slotrectangular spouted bed, 7th Int. Conf. Discret. Elem. Methods 2016, pp. 697–707, https://doi.org/10.1007/978-981-10-1926-5. [57] F. Marchelli, C. Moliner, B. Bosio, E. Arato, A CFD-DEM sensitivity analysis: the case of a pseudo-2D spouted bed, Powder Technol. 353 (2019) 409–425, https://doi.org/10. 1016/j.powtec.2019.05.035. [58] C. Moliner, F. Marchelli, L. Ong, A. Martinez-Felipe, D. Van der A, E. Arato, Sensitivity analysis and validation of a Two Fluid Method (TFM) model for a spouted bed, Chem. Eng. Sci. 207 (2019) 39–53, https://doi.org/10.1016/j.ces.2019.06.008. [59] L.W. Rong, K.J. Dong, A.B. Yu, Lattice-Boltzmann simulation of fluid flow through packed beds of uniform spheres: effect of porosity, Chem. Eng. Sci. 99 (2013) 44–58, https://doi.org/10.1016/J.CES.2013.05.036. [60] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Géotechnique. 29 (1979) 47–65, https://doi.org/10.1680/geot.1979.29.1.47. [61] Y. Tsuji, T. Tanaka, T. Ishida, Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe, Powder Technol. 71 (1992) 239–250, https://doi.org/10.1016/0032-5910(92)88030-L. [62] Y. Tsuji, T. Kawaguchi, T. Tanaka, Discrete particle simulation of two-dimensional fluidized bed, Powder Technol. 77 (1993) 79–87, https://doi.org/10.1016/00325910(93)85010-7. [63] C. Moliner, F. Marchelli, N. Spanachi, A. Martinez-Felipe, B. Bosio, E. Arato, CFD simulation of a spouted bed: comparison between the discrete element method (DEM) and the two fluid model (TFM), Chem. Eng. J. 377 (2019) 120466, https://doi.org/10. 1016/j.cej.2018.11.164. [64] S. Golshan, B. Esgandari, R. Zarghami, CFD-DEM and TFM simulations of spouted bed, Chem. Eng. Trans. 57 (2017) 1249–1254, https://doi.org/10.3303/CET1757209. [65] A. Stroh, F. Alobaid, M.T. Hasenzahl, J. Hilz, J. Strohle, B. Epple, Comparison of three different CFD methods for dense fluidized beds and validation by a cold flow experiment, Particuology. 29 (2016) 34–47, https://doi.org/10.1016/j.partic.2015.09.010. [66] X.-L. Zhao, S.-Q. Li, G.-Q. Liu, Q. Yao, J.-S. Marshall, DEM simulation of the particle dynamics in two-dimensional spouted beds, Powder Technol. 184 (2008) 205–213, https://doi.org/10.1016/J.POWTEC.2007.11.044. [67] G.-Q. Liu, S.-Q. Li, X.-L. Zhao, Q. Yao, Experimental studies of particle flow dynamics in a two-dimensional spouted bed, Chem. Eng. Sci. 63 (2008) 1131–1141, https:// doi.org/10.1016/J.CES.2007.11.013. [68] M. Olazar, M.J. San José, M.A. Izquierdo, S. Alvarez, J. Bilbao, Local bed voidage in spouted beds, Ind. Eng. Chem. Res. 40 (2001) 427–433, https://doi.org/10.1021/ ie0003741. [69] M. Olazar, M.J. San José, M.A. Izquierdo, A. Ortiz De Salazar, J. Bilbao, Effect of operating conditions on solid velocity in the spout, annulus and fountain of spouted beds, Chem. Eng. Sci. 56 (2001) 3585–3594, https://doi.org/10.1016/S0009-2509 (01)00022-7. [70] D. Geldart, Types of gas fluidization, Powder Technol. 7 (1973) 285–292, https://doi. org/10.1016/0032-5910(73)80037-3. [71] ANSYS, Chapter 1: Basic fluid flow, ANSYS FLUENT Theory Guid, ANSYS, 2018. [72] ANSYS, Chapter 16: Discrete phase, ANSYS FLUENT Theory Guid, ANSYS, 2018. [73] ANSYS, Chapter 17: Multiphase flows, ANSYS FLUENT Theory Guid, ANSYS, Inc, 2019. [74] J. Xie, W. Zhong, B. Jin, LES-Lagrangian modelling on gasification of combustible solid waste in a spouted bed, Can. J. Chem. Eng. 92 (2014) 1325–1333, https://doi. org/10.1002/cjce.21941. [75] C. Hu, K. Luo, S. Wang, L. Junjie, J. Fan, The effects of collisional parameters on the hydrodynamics and heat transfer in spouted bed: a CFD-DEM study, Powder Technol. 353 (2019) 132–144, https://doi.org/10.1016/J.POWTEC.2019.05.020.

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058

18

F. Marchelli et al. / Powder Technology xxx (2019) xxx

[76] S.I. Rubinow, J.B. Keller, The transverse force on a spinning sphere moving in a viscous fluid, J. Fluid Mech. 11 (1961) 447, https://doi.org/10.1017/ S0022112061000640. [77] R.-J. Liu, R. Xiao, M. Ye, Z. Liu, Analysis of particle rotation in fluidized bed by use of discrete particle model, Adv. Powder Technol. 29 (2018) 1655–1663, https://doi. org/10.1016/J.APT.2018.03.032. [78] Q. Zhou, L.-S. Fan, Direct numerical simulation of low-Reynolds-number flow past arrays of rotating spheres, J. Fluid Mech. 765 (2015) 396–423, https://doi.org/10. 1017/jfm.2014.732. [79] Y. Wang, L. Zhou, Q. Yang, Hydro-mechanical analysis of calcareous sand with a new shape-dependent fluid-particle drag model integrated into CFD-DEM coupling program, Powder Technol. 344 (2019) 108–120, https://doi.org/10.1016/J.POWTEC. 2018.12.008. [80] G. Bagheri, C. Bonadonna, On the drag of freely falling non-spherical particles, Powder Technol. 301 (2016) 526–544, https://doi.org/10.1016/J.POWTEC.2016.06.015. [81] E. Loth, Drag of non-spherical solid particles of regular and irregular shape, Powder Technol. 182 (2008) 342–353, https://doi.org/10.1016/J.POWTEC.2007.06.001. [82] W. Zhong, A. Yu, X. Liu, Z. Tong, H. Zhang, DEM/CFD-DEM modelling of non-spherical particulate systems: theoretical developments and applications, Powder Technol. 302 (2016) 108–152, https://doi.org/10.1016/j.powtec.2016.07.010.

[83] B. Ren, W. Zhong, Y. Chen, X. Chen, B. Jin, Z. Yuan, Y. Lu, CFD-DEM simulation of spouting of corn-shaped particles, Particuology. 10 (2012) 562–572, https://doi. org/10.1016/j.partic.2012.03.011. [84] X. Liu, J. Gan, W. Zhong, A. Yu, Particle shape effects on dynamic behaviors in a spouted bed: CFD-DEM study, Powder Technol. (2019) https://doi.org/10.1016/J. POWTEC.2019.07.099. [85] D. Gidaspow, R. Bezburuah, J. Ding, Hydrodynamics of circulating fluidized beds, kinetic theory approach, Fluid. VII - Proc. 7th Eng. Found. Conf. Fluid 1992, pp. 75–82. [86] S. Tenneti, R. Garg, S. Subramaniam, Drag law for monodisperse gas–solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres, Int. J. Multiphase Flow 37 (2011) 1072–1092, https://doi.org/10.1016/J. IJMULTIPHASEFLOW.2011.05.010. [87] D. Liu, C. Bu, X. Chen, Development and test of CFD–DEM model for complex geometry: a coupling algorithm for Fluent and DEM, Comput. Chem. Eng. 58 (2013) 260–268, https://doi.org/10.1016/j.compchemeng.2013.07.006. [88] ANSYS, Chapter 21: Solver theory, ANSYS FLUENT Theory Guid. ANSYS, 2018. [89] F. Marchelli, D. Bove, C. Moliner, B. Bosio, E. Arato, Discrete element method for the prediction of the onset velocity in a spouted bed, Powder Technol. 321 (2017) 119–131, https://doi.org/10.1016/j.powtec.2017.08.032.

Please cite this article as: F. Marchelli, Q. Hou, B. Bosio, et al., Comparison of different drag models in CFD-DEM simulations of spouted beds, Powder Technol., https://doi.org/10.1016/j.powtec.2019.10.058