Powder Technology 264 (2014) 343–364 A comparative study on the influence of the gas flow rate on the hydrodynamics of a gas-solid spouted fluidized bed using Euler-Euler and Euler-Lagrange / DEM models Naser Almohammeda,∗, Falah Alobaidb , Michael Breuera , Bernd Eppleb a b Professur für Strömungsmechanik, Helmut-Schmidt-Universität Hamburg, Holstenhofweg 85, D-22043 Hamburg, Germany Institut für Energiesysteme und Energietechnik, Technische Universität Darmstadt, Otto-Berndt-Straße 2, D-64287 Darmstadt, Germany Abstract In this work, different modeling approaches are used to study the hydrodynamics of gas-solid flows in a three-dimensional, lab-scale spouted fluidized bed. In particular, the simulation results obtained by the two-fluid model are compared with those of coupled CFD/DEM simulations. To explore the effects of the gas mass flow rate on the ability of the used simulation techniques to predict the flow behavior in the simulated test case, two different fluidization conditions are considered while maintaining the same computational set-up. The spouted fluidized bed has been evaluated by means of high-speed imaging for validating the simulation results. A comparative study of the two modeling approaches with experimental observations for bubble size in terms of the representative bubble diameter and the bed height has been carried out under both operating conditions. To identify significant simulation settings for the two-fluid model, a parameter study is carried out. The investigated input parameters include the gas-particle drag models, the granular temperature approach, the solid-phase wall boundary conditions in terms of specularity coefficient and the particle-particle restitution coefficient. At a gas mass flow rate of 0.005 kg/s, the predicted characteristics of the bubble formation and the bed expansion using both techniques are in very good agreement with experimental observations. Using the same validated settings, the two-fluid model shows some notable discrepancies when the gas mass flow rate is increased to 0.006 kg/s, while the DEM is more successful than the two-fluid model in reproducing realistic flow patterns. Although both techniques predict the right fluidization regimes and trends in bubble sizes and bed expansion, their results deviate significantly from the experimental data during the final stage of the bubble formation. The most important reasons for the differences in predicting the bed dynamics along with the advantages and limitations of the two approaches are discussed. Keywords: Spouted Fluidized Bed, Euler-Euler Model, Euler-Lagrange Model, Two-Fluid Model, Discrete Element Method 1. Introduction A gas-solid fluidized bed is typically a packed bed of particulate solids in a container with perforated base, through which an upward flow of gas is passed. Fluidized bed reactors have been widely applied to various industrial processes including the coating, the granulation, the combustion and the gasification of coal, biomass and solid waste. Here, the corresponding multiphase flows, which are of substantial interest in the recent academic and industrial research, play a significant role. The modeling, the characterization and the advanced understanding of the hydrodynamics of these complex flows are still challenging because their behavior is still unknown and difficult to be predicted in many applications. Therefore, the basic objective of modeling the hydrodynamics of the particulate flow behavior in a fluidized bed is to obtain a realistic flow structure. Experiments have provided substantial insights into the hydrodynamic behavior of particulate flows in lab-scale fluidized beds, but they are difficult, expensive, and limited to this scale of fluidized beds. ∗ Corresponding author. Tel.: +49 40 65412291; Fax: +49 40 65413781 Email address: naser.almohammed@hsu-hh.de (Naser Almohammed) Powder Technology 264 (2014) 343–364 On the other hand, the analysis of the theoretical approaches proposed in the past was restricted to simplified models. In recent years, together with rapid development of high-performance computers and numerical algorithms, computational fluid dynamics (CFD) has proved to be a powerful tool to bridge the gaps between analytical solutions and experiments. Although the applications of CFD in dispersed multiphase flows are still limited when compared to single-phase flows due to their complex flow behavior, advances in computational fluid dynamics, however, have provided the basis for further insight into the hydrodynamics and related characteristics of multiphase flow systems. In the last decade, significant research efforts have been made to develop detailed models to investigate the particulate flow in gas-solid fluidized beds. A review on models used to study gas-fluidized beds can be found, among other publications, in Van der Hoef et al. (2008), Deen et al. (2007) and Zhu et al. (2008). Broadly speaking, two hydrodynamic models are distinguished: the Euler-Euler approach, or commonly called two-fluid (continuum) model, and the Euler-Lagrange approach (discrete particle model). These techniques have been commonly adopted in the CFD modeling of the gas-solid flow in fluidized beds. Both models consider the gas phase as a continuum, and the main difference is how to deal with the dispersed phase. The basic idea of the two-fluid model (TFM) is that both phases are treated as fully interpenetrating continua, and the challenge facing the modeling of the motion of the solid phase is how to model the kinetic and collisional transport of the particles (gas-particle, particle-particle and particle-wall interactions) while representing the granular phase as an interpenetrating continuum. Therefore, many assumptions have been made regarding the gas-particle drag correlations and the rheology of the solid phase. Euler-Lagrange models on the other hand do not require additional closure equations for the particulate phase since the motion of each individual element (particle or droplet) is computed based on first principles. However, the modeling of dense gas-solid flows based on the discrete particle model (DPM) requires an accurate description of the particle-particle and particle-wall collisions. These collisions can be stochastically modeled or deterministically detected. In dense gas-solid flows, stochastic collision models produce, however, unrealistic results such as solid volume fractions greater than unity in the control volume (Götz, 2006). Therefore, these models are mainly used for the dilute gas-solid flows. For the discrete particle simulations with deterministic collision detection, two models are widely applied: the hard-sphere model and the soft-sphere model. In the hard-sphere model, single binary collisions are considered as instantaneous processes. In the soft-sphere model, also known as discrete element method (DEM), the particles can overlap each other or penetrate into the wall. Depending on the penetration depth, a contact force changing the motions of particles is determined (Cundall and Strack, 1979; Tsuji et al., 1992, 1993). Although the deterministic collision models show much better agreement with the experimental data than stochastic collision models, they are typically impractical for simulating industrial systems owing to their high computational demands. The stochastic collision models represent today the only possibility to numerically calculate the practically relevant gas-solid flows in large-scale facilities due to their high efficiency. For both simulation techniques adopted in this work, previous studies relevant to this topic are briefly summarized below according to the modeling approach and are followed by the original contribution of this work. 1.1. Two-Fluid Model In the framework of the two-fluid model, many studies have been conducted to predict the hydrodynamic behavior of gas-solid flows in cold-flow spouted fluidized beds. Huilin et al. (2004) combined a hydrodynamic two-fluid model with a kinetic-frictional constitutive model of solids to predict the hydrodynamics of the gas-solids flow in a spouted bed. They concluded that the formation of the jet zone, the fountain zone and the annulus zone depends significantly on the variation of the frictional stresses. In addition, changing the values of the particle-particle restitution coefficient results in different flow structures in the jet zone. Du et al. (2006a,b) investigated the impact of several gas-solid drag models, the frictional stress, the maximum packing limit and the restitution coefficient on the simulation of spouted fluidized beds based on the two-fluid model. They found that the hydrodynamics of spouted beds depends strongly on the restitution coefficient since changing its value also affects the granular temperature, which represents the fluctuating energy of the solids. Duarte et al. (2009) used the Euler-Euler model to simulate the hydrodynamic behavior in conical and conicalcylindrical spouted beds. Gryczka et al. (2009) modeled the hydrodynamic behavior of a two-dimensional gas-solid flow in a prismatic apparatus adapting the two-fluid model. In their study, various gas-particle drag correlations, two granular temperature approaches and different assumptions for the value of the particle-particle restitution coefficient were tested. Furthermore, the characteristic gas fluctuations over the entire bed were recorded by a high frequency pressure detector and compared with calculated fluctuations. Furthermore, a Fourier analysis on the measured and simulated pressure 2 Powder Technology 264 (2014) 343–364 frequency spectra were contrasted. Lan et al. (2012) investigated the influence of solid-phase wall boundary conditions in terms of specularity coefficient and particle-wall restitution coefficient on CFD simulations of spouted beds. Their simulated results showed that the solid-phase wall boundary conditions play an important role in the CFD simulation of spouted beds, and the specularity coefficient has a significant effect on the fountain height, the wall shear stress, the pressure drop and the velocity distribution of the solids. Furthermore, they observed that the particle-wall restitution coefficient plays only a minor role in the modeling of the gas-solid flow in spouted beds. Pei et al. (2011) applied the two-fluid model with a zero particle viscosity to investigate the effects of inter-phase drag force on the jet characteristics in a two-dimensional spouted fluidized bed. In their work, the obtained simulation results using different drag force models were compared with the experiment. 1.2. Discrete Element Method The hydrodynamics of the gas-solid flow in the fluidized bed adopting the CFD/DEM model has been investigated in several studies: Link (2006) developed and validated a discrete particle model of a spouted-fluid bed granulator. Different closure laws of the gas-particle drag were tested, and the simulation results concerning gas phase pressure fluctuations over the entire bed were compared with the experimental results in a pseudo-2D spouted bed. Although the simulations showed satisfactory agreement with the measurements, the CFD/DEM calculation is computationally very expensive. Götz (2006) and Tsuji and Tanaka (2008) parallelized the CFD/DEM simulation by decomposing the computational domain into several sub-blocks. It was concluded that an increase in the number of processors results in an almost linear increase of the performance. When exceeding the critical number of the decomposed sub-blocks, the performance is no longer proportional to the number of processors due to the fact that the communication time for the data exchange between the processors increases massively. Alobaid and Epple (2013) investigated the hydrodynamic behavior of the dense gas-solid flow in a lab-scale spouted fluidized bed employing an in-house CFD/DEM code. Here, the fluid-particle interaction was performed using a new procedure called the offset method in which the grid cells undergo several numbers of spatial displacements leading to improved calculation of variables (e.g., the volume fraction) concerning the momentum exchange between the two phases. Alobaid et al. (2013) introduced a new procedure, which allows the variation of the fluid grid resolution independent of the particle size. To achieve this, they introduced an additional grid, the so-called particle grid, in which the physical values of the solid phase is computed. A validation study was carried out to assess the obtained numerical results during 0.5 s operation of a Plexiglas spouted-fluidized bed. A good qualitative correlation of the solid distribution in the bed and the acceptable quantitative agreement of the pressure drops at different positions along the bed height were obtained. Besides the above mentioned works, Gera et al. (1998) used the DEM and TFM to study the hydrodynamics of a two-dimensional fluidized bed at one fluidization velocity. They found that the predictions of both techniques are in good qualitative agreement with the experiment. Furthermore, they discussed the quantitative differences of both approaches for a single bubble rising in a two-dimensional fluidized bed. Goldschmidt et al. (2004) incorporated the two-fluid model and the hard-sphere model as an Euler-Lagrange approach to investigate the hydrodynamic behavior of the dense particulate flow in the spouted fluidized beds. They validated and compared a three-dimensional hard-sphere model as well as a continuum model based on the kinetic theory of granular flow (KTGF). Their study indicated that the limited agreement of the hydrodynamic model with discrete particle simulations and experiments occurs because the rotation of the particles is not taken into account in the two-fluid model. In addition, a better agreement between the results obtained employing the hard-sphere model and the experiments was achieved by an effective restitution coefficient. Chiesa et al. (2005) compared the prediction of bubble formation in a two-dimensional gas fluidized bed by both the Euler-Lagrange approach with a deterministic hard-sphere model and the Euler-Euler approach. Their numerical predictions were furthermore compared to experimental observations of the bubble formation in the bed. They concluded that the results from both approaches are from a macroscopic point of view satisfactory and the numerical simulations show a rather similar pattern to what was observed in the experiment. The results obtained by the EulerLagrange approach showed a much better agreement with the experimental results than those obtained by the continuous approach. However, the CPU time necessary to perform the numerical simulation with the Euler-Lagrange model was four orders of magnitude higher than the time required to perform an Euler-Euler simulation. 3 Powder Technology 264 (2014) 343–364 1.3. The objectives of the paper The above mentioned studies suggest that both models can predict the complex gas-solid flows in different types of spouted beds, but they did not investigate the capability of both techniques to reproduce the flow patterns when the flow rate of the gas passing through the bed is changed using the same computational set-up. To the best of our knowledge, there has been no previous study dealing with the comparison between the two-fluid model and the discrete element method against experiments at different fluidization conditions. Therefore, the main goal of the present study is to investigate the hydrodynamic behavior of three-dimensional gas-solid flows in a spouted fluidized bed model consisting of monodisperse particles at two different gas mass flow rates using the two-fluid model and the discrete element method. For this purpose, new experiments are carried out at gas mass flow rates of 0.005 kg/s and 0.006 kg/s, respectively. Starting with a gas mass flow rate of 0.005 kg/s, the influence of the different simulation parameters of the two-fluid model on the overall bed behavior is tested. Then, the best simulation results of both techniques are compared against each other and validated with the experiments by analyzing the bubble size in terms of the equivalent bubble diameter and the bed height during 500 ms. Finally, the gas mass flow rate is increased to 0.006 kg/s and the same validated computational set-up is used to simulate the same fluidized bed, and the obtained numerical results of both techniques are compared again with the measurements. This paper is organized in the following manner. First, the two-fluid model and the discrete element method adopted in this work are briefly explained. Then, the experiment and numerical set-up of both techniques are presented. A detailed comparison of the obtained results using both models with the experiments is reported. Finally, conclusions concerning the accuracy and efficiency of these models are drawn. 2. Mathematical and Numerical Methods 2.1. Euler-Euler Model (Two-Fluid Model) In the Euler-Euler approach, two phases are treated mathematically as fully interpenetrating continua by introducing the concepts of the solid-phase viscosity and the solid-phase pressure. 2.1.1. Governing equations This section presents the governing equations of the two-fluid model available in FLUENT 14.0 (ANSYSr , 2012) and used in this comparative analysis. In this work, two phases are considered: a gas phase (air) and a granular solid phase (glass beads). The particles of the granular solid phase are assumed to be spherical and monodisperse. The mass conservation equation for the phase q is written as: n X ∂ (ṁ pq − ṁqp ) . (αq ρq ) + ∇ · (αq ρq~uq ) = ∂t p=1 (1) In the above equation, αq stands for the phasic volume fraction, which represents the space occupied by each phase, ρq and ~uq are the density and velocity vector of the phase q, respectively. The symbol ṁ pq characterizes the inter-phase mass transfer from the pth to the qth phase, and ṁqp denotes the inter-phase mass transfer from phase q to phase p. The total volume fraction, which is the sum of all the volume fractions, in each computational cell must be equal to n P αq = 1. In this study, no inter-phase mass transfer takes place between both phases. Hence, the inter-phase unity q=1 mass transfer from the solid phase (particles) s to the continuous phase (gas) g and vice versa are set equal to zero (ṁ sg = ṁgs = 0). The momentum conservation equation for the phase q has the following form (Gryczka et al., 2009): ∂ (αq ρq~uq ) + ∇ · (αq ρq~uq~uq ) = − αq ∇pq + ∇ · τq ∂t   + αq ρq F~ex,q + F~lift,q + F~vm,q n  X  ~ pq + ṁ pq~u pq − ṁqp~uqp , R + p=1 4 (2) Powder Technology 264 (2014) 343–364 where F~ex,q , F~lift,q and F~vm,q are the external body force, the lift force and the virtual mass force, respectively. In most cases, the lift force is insignificant compared to the drag force, so there is no reason to include this extra term. In this work, the lift force is not considered because the inclusion of this force is not appropriate for closely packed particles (ANSYSr , 2012). Furthermore, the Basset force, the virtual mass force and the pressure gradient force are assumed to be negligible compared to the drag fore when the solid phase density is much higher than the gas phase density ρ p /ρg ≫ 1 (Hjelmfelt and Mockros, 1966; Armenio and Fiorotto, 2001). It is well known that, in a fluidized bed, the relative acceleration or deceleration rate of a particle is low. Hence, the characteristic time scale of these changes is much smaller than the Stokes relaxation time. Therefore, only the influence of the interaction force (or drag force) ~ pq and the gravity force are considered in Eq. (2). The symbols ~u pq and ~uqp stand for the inter-phase between the phases R velocities which are defined as ~u p or ~uq depending on the inter-phase mass transfer. However, since no inter-phase mass transfer is taken into account here, the corresponding terms in Eq. (2) cancel out. τq stands for the stress-strain tensor of the qth phase and is defined as (Gryczka et al., 2009): !   2 τq = αq µq ∇~uq + ∇~uTq + αq λq − µq ∇ · ~uq I , (3) 3 where µq and λq express the shear and bulk viscosity of the phase q, respectively. I is the unit tensor. The momentum ~ pq , which is subject conservation given by Eq. (2) must be closed by appropriate expressions for the inter-phase force R ~ ~ ~ to the conditions that R pq = −Rqp and Rqq = 0. In the two-fluid model, a simple interaction term of the following ~ pq = K pq (~u p − ~uq ), where K pq (= Kqp ) is the inter-phase momentum exchange coefficient. The drag force form is used R acting on a single particle in a gas-particle system can be expressed by the product of a gas-solid momentum exchange ~ gs = β (~ug − ~u s ). coefficient β and the relative (or slip) velocity between the gas and solid phases R 2.1.2. Drag models The correlations for the momentum transfer coefficient β are usually obtained from the pressure drop measurements in fixed, fluidized or settling beds. Numerous correlations (drag models) for calculating the momentum exchange coefficient of gas-solid systems have been reported in the literature. The momentum exchange coefficient between the gas phase g and the solid phase s can be written as: β= αsρs f . τs (4) In the above equation, f denotes the drag function, and τ s is the so-called relaxation time of particles given by: τs = ρ s d2s , 18µg (5) where d s denotes the diameter of particles in the solid phase and µg is the gas viscosity. In this work, the gas-solid momentum exchange coefficient is calculated by applying three different drag models among which only the drag function is different. These drag models are: • Syamlal and O’Brien (1989): In this model, the drag function is expressed as: f = C D Re s αg 24v2r,s , (6) where C D stands for the drag coefficient and is derived by Dalla Valle (1948): 2   4.8   C D = 0.63 + p  , Re s /vr,s and vr,s is the terminal velocity correlation for the solid phase defined as: p   vr,s = 0.5 A − 0.06 Re s + (0.06 Re s )2 + 0.12 Re s (2B − A) + A2 5 (7) (8) Powder Technology 264 (2014) 343–364 1.28 for α ≤ 0.85 and B = α2.65 for α > 0.85. The symbol Re stands for the solids with A = α4.14 g g s g , and B = 0.8 αg g Reynolds number defined as: Re s = ρg d s ~u s − ~ug µg . (9) This model is based on the measurements of the terminal velocities of particles in fluidized or settling beds. Thus, the resulting gas-solid exchange coefficient can be written as: ! 3 α s αg ρg Re s ~u s − ~ug . β = CD 2 (10) 4 vr,s d s vr,s • Wen and Yu (1966): In this model, the gas-solid exchange coefficient is of the following form: 3 α s αg ρg ~u s − ~ug −2.65 αg β = CD 4 ds (11) here C D is calculated based on the following equation: CD = i 24 h 1 + 0.15(αg Re s )0.687 αg Re s (12) This model is appropriate for dilute systems. • Gidaspow et al. (1992): This model is a combination of the Wen and Yu (1966) model and the Ergun (1952) equation. For a solid volume fraction αg > 0.8 (i.e., very low particle concentrations), the gas-solid exchange coefficient β is calculated from the previous model of Wen and Yu (1966). For αg ≤ 0.8, the gas-solid exchange coefficient β is calculated by Ergun’s equation based on the experimental data of the pressure drop over fixed, dense beds of monodisperse particles:   α s µg 1 − αg ρg α s ~u s − ~ug . (13) + 1.75 β = 150 2 ds αg d s This model is recommended for dense fluidized beds. 2.1.3. Kinetic theory of granular flow In dense gas-solid flows, the particles are moved by the kinetic transport due to the interaction with the gas as well as a collisional transport due to the particle-particle collisions. In the two-fluid model, the modeling of the particle motion is based on the kinetic theory for dense gases. The kinetic theory of granular flow (KTGF) is basically an extension of the classical kinetic theory of gases to dense particulate flows. In this theory, non-ideal particle-particle collisions and the gas-particle drag are taken into account. Analogous to the kinetic theory of gases, particles are allowed to travel freely ~ s and a and to collide with neighboring particles. The particle velocity can be decomposed into a mean local velocity U ~ ~ ~ superimposed fluctuating component C s representing the vibration of the particle, so that ~u s = U s + C s . The granular temperature is introduced to represent the solids fluctuating energy. Analogous to the thermodynamic temperature of the gas, the granular temperature Θ s is associated with the random velocity fluctuations of the solid particles: Θs = 1 D~ ~ E Cs · Cs . 3 (14) The transport equation of the granular temperature Θ s derived from the kinetic theory takes the form: ! 3 ∂ (α s ρ s Θ s ) + ∇ · (α s ρ s~u s Θ s ) = (−p s I + τ s ) : ∇~u s 2 ∂t
+ ∇ · ΓΘs ∇Θ s − γΘs + φls , 6 (15) Powder Technology 264 (2014) 343–364 where (−p s I + τ s ) : ∇~u s stands for the generation of energy by the solid stress tensor, ΓΘs ∇Θ s denotes the diffusion of energy, γΘs is the collisional dissipation of energy and φls represents the energy exchange between the lth fluid or solid phase and the sth solid phase. According to Gidaspow et al. (1992), the diffusion coefficient for granular energy ΓΘs is expressed as: √ #2 " 6 150ρ s d s Θ s π ΓΘs = 1 + α s g0,ss (1 + e ss ) 384 (1 + e ss ) g0,ss 5 (16) r Θs 2 + 2ρ s α s d s (1 + e ss )g0,ss . π The rate of energy dissipation within the sth solids phase due to inelastic collisions is given by the expression derived by Lun et al. (1984):   12 1 − e2ss g0,ss ρ s α2s Θ3/2 (17) γΘ s = √ s , ds π and φgs = −3 β Θ s . (18) In the above equations, e ss denotes the restitution coefficient of particle-particle collisions, and g0,ss stands for the radial distribution function, which is a correction factor modifying the probability of collisions between grains within the granular phase when it becomes dense. g0,ss is defined as:  !1/3 −1   α s g0,ss = 1 − (19)  . α s,max The function g0,ss governs the transition from the “compressible” condition with α s < α s,max , where the spacing between the solid particles can continue to decrease, to the “incompressible” condition with α = α s,max , where no further decrease in the spacing can occur. The granular temperature is used to calculate the solids pressure and viscosity and the particle wall boundary conditions explained next. 2.1.4. Solids pressure and viscosity In the two-fluid model, the solids pressure is composed of a kinetic term and a second term due to particle collisions. p s = α s ρ s Θ s + 2ρ s (1 + e ss ) α2s g0,ss Θ s (20) To calculate the stress-strain tensor of the solid phase given by Eq. (3), the shear and bulk viscosity (i.e., µ s and λ s ) are still unknown. The granular viscosity is a summation of three viscosity contributions µ s,col , µ s,kin and µ s,fr , which stand for the collisional and kinetic parts and the optional frictional part. They are added to give the solids shear viscosity: µ s = µ s,col + µ s,kin + µ s, f r (21) with: • The collisional viscosity by Gidaspow et al. (1992): r Θs 4 2 µ s,col = α s ρ s d s g0,ss (1 + e ss ) 5 π (22) • The kinetic viscosity by Gidaspow et al. (1992): √ " #2 4 10ρ s d s Θ s π 1 + g0,ss α s (1 + e ss ) α s µ s,kin = 96α s (1 + e ss ) g0,ss 5 7 (23) Powder Technology 264 (2014) 343–364 • The frictional viscosity by Schaeffer (1987): µ s, f r = p s sin φ √ 2 I2D (24) where φ is the angle of internal friction, and I2D is the second invariant of the deviatoric stress tensor. • The solids bulk viscosity is expressed by Lun et al. (1984) as: r 4 Θs . λ s = α s ρ s d s g0,ss (1 + e ss ) 3 π (25) 2.1.5. Particle wall boundary conditions At the wall, the tangential component of the velocity of particles and the granular temperature are defined applying the Johnson and Jackson (1987) wall boundary conditions. p π√ αs τs = − ρ p g0,ss Θ s u s,w , 3ϕ (26) 6 α s,max where ϕ stands for the specularity coefficient. The slip velocity between particles and the wall u s,w can be obtained by assuming that the tangential force exerted on the boundary and the force due to the particles shear stress close to the wall are equal (Lan et al., 2012): u s,w = − √ ∂u s,w . 3Θ s πϕρ s α s g0,ss ∂n 6µ s α s,max (27) In a similar manner, the granular temperature at the wall is also obtained by assuming that the granular temperature flux at the wall due to the energy dissipation and the generation of granular energy due to the slip at the wall region are equal (Lan et al., 2012): √ 3Θπϕρ s α s u2s,slip g0,ss Θ3/2 s k s Θ s ∂Θw + (28) Θw = − γw ∂n 6α p,max γw with √ 3π(1 − e2w )ρ s α p g0,ss Θ3/2 s , γw = 4α s,max (29) where ew is the particle-wall restitution coefficient and k s stands for the turbulent kinetic energy of the solids. The specularity coefficient ϕ is an empirical parameter used to quantify the nature of particle-wall collisions. In other words, it characterizes the friction between the particles and the wall. The values of the specularity coefficient are in the interval 0 ≤ ϕ ≤ 1. A value of zero implies a specular collision (free-slip conditions) and unity corresponds to diffuse collision (no-slip conditions). Thus, the nature of the particle-wall collisions can be expected to influence the shear stress and the fluctuating energy flux at the boundary. 2.1.6. Turbulence model In this work, a the k − ǫ dispersed turbulence model is used (Gryczka et al., 2009; ANSYSr , 2012). Here, the turbulent viscosity is defined as: µt,g = ρgCµ kg2 ǫg , (30) where ǫg stands for the dissipation rate of the turbulent kinetic energy of the gas phase kg , and the constant Cµ = 0.09. The transport equation of the turbulent kinetic energy of the gas phase is expressed as: ! µt,g ∂ (αg ρg kg ) + ∇ · (αg ρg~ug kg ) =∇ · αg ∇kg + αgGk,g ∂t σk (31) − αg ρg ǫg + αg ρg Πkg , 8 Powder Technology 264 (2014) 343–364 and the transport equation of the dissipation rate of the turbulent kinetic energy of the gas phase is expressed as: ! µt,g ǫg ∂ (αg ρg ǫg ) + ∇ · (αg ρg~ug ǫ) =∇ · αg ∇ǫg + αg C1ǫ Gk,g ∂t σǫ kg − αg ǫg2 kg (32) C2ǫ ρg + αg ρg Πǫg . In the above two equations, Gk,g is the production of turbulent kinetic energy, and the other terms have the same meaning as in the single-phase k − ǫ model. Πkg and Πǫg represent the influence of the dispersed phases on the continuous phase, and the relationship between them is modeled according to Elghobashi and Abou-Arab (1983) as: Πǫg = C3ǫ ǫg Πk kg g (33) with C1ǫ = 1.44, C2ǫ = 1.92, C3ǫ = 1.2, σk = 1 and σǫ = 1.3. The turbulent drag term in Eq. (2) is modeled as follows (Gryczka et al., 2009; ANSYSr , 2012): ~ pq = K pq (~u p − ~uq ) = (U ~p − U ~ q ) − K pq~udr , R (34) where ~udr stands for the drift velocity resulting from the turbulent fluctuations in the volume fraction. When multiplied by the exchange coefficient K pq , it serves as a correction to the momentum exchange term for turbulent flows. The drift velocity is given by: ! Dq Dp ~udr = − (35) ∇α p − ∇αq , σ pq α p σ pq αq with a dispersion Prandtl number σ sg = 0.75. Here, it is assumed that the diffusivities are equal, i.e., D p = Dq = Dt,pq ). The diffusivity is modeled by: 1 Dt,pq = k pq τt,pq , 3 (36) where τt,pq denotes the Lagrangian integral time scale calculated along particle trajectories. 2.2. Euler-Lagrange Model (Discrete Element Method) Beside the Eulerian-Eulerian approach, the Eulerian-Lagrangian approach with deterministic collision models are a growing alternative model for the numerical representation of gas-solid flows. The DEM model offers detailed information on the hydrodynamic behavior of granular flows. Here, the determination of a single particle trajectory and the angular displacement is achieved by establishing balances of forces and moments of force on each particle with the aid of Newton’s equations of motion. 2.2.1. Kinetics of the particles The translational and rotational motions of a particle (solid phase), denoted here by p, are calculated by the integration of Newton’s equations of motion: d~u p X F~ p,k = dt k d~ ωp X ~ p,k Ip M = dt k mp (37) ~ p are the particle translational and angular velocities. The symbols m p and I p are the particle mass and where ~u p and ω ~ p,k are the force and the moment of force. The index k denotes the the moment of inertia, respectively. F~ p,k and M 9 Powder Technology 264 (2014) 343–364 number of acting forces or moments on the particle. In this study, the considered fluid force acting on a particle F~ p is given by n F~ p = F~gra + F~ pre + F~dra + F~ sa f + F~mag + F~con + F~ t | {z con } (38) Contact forces In the above equation, F~gra is the gravitational force, F~ pre the pressure gradient (or buoyancy) force, F~dra the drag force, F~ sa f the Saffman lift force, F~mag the Magnus force and F~con the forces resulting from the interactions between the particles themselves and the wall. In this section, only the calculation of the contact force is explained, and the remaining applied forces are listed in Appendix A. Note that the inter-particle forces such as the van der Waals force and the electrostatic force are not considered in the present DEM simulations because the used particles are relatively large (d p = 2.5 mm). Furthermore, the particles are assumed to be dray and electrostatically neutral. Hence, the attractive forces with physical connection (e.g. liquid bridge and sintering bridge) are also omitted. 2.2.2. Particle–Particle interaction The considered contact between the two spherical particles i and j with the position vectors ~ri and ~r j , the radii ~ i and ω ~ j is shown in Figure 1. The ri and r j , the masses mi and m j , the translational and angular velocities ~ui , ~u j , ω determination of the contact force due to the particle-wall collision can be treated analogously, however, under the assumptions that the wall has an infinite radius and the value of the translational and angular velocities are zero. ωj ωi u Particle j Particle i u u Wall Particle i i,w ωi (a) i,w (b) Figure 1: DEM: Contact forces due to (a) particle-particle collision and (b) particle-wall collision. In the DEM model, the particles can overlap each other or penetrate into the wall. Depending on the penetration depth, the contact force F~con can be modeled using a spring-damper-slider system. The resulting contact force on the particle i is obtained by summing up all the normal and tangential contact forces that act simultaneously between the particle i and other particles as well as the walls: F~con,i = N  X k=1 k,i  n t . F~con,ik + F~con,ik (39) The index N represents the number of collision partners (i.e., particles and walls that are in contact with the particle i). n (Hertz, 1882) Based on the Voigt-Kelvin model, the normal contact force results from the sum of the elastic force F~ela n and the damping force F~dam and is described by the differential equation of the spring-damper system: n F~con,i j = mi j d~unij dt
= −kn δn 3/2 ~ni j − ηn~unij , | {z } |{z} n F~ela (40) n F~dam 10 Powder Technology 264 (2014) 343–364 with the normal unit vector ~ni j , which points in the direction from the center of the particle j to the center of the particle i as depicted in Figure 1, defined as: ~ri − ~r j ~ni j = ~ri − ~r j . (41) In Eq. (40), kn is the normal stiffness constant, ηn is the normal coefficient of damping, mi j is the reduced mass (effective mass) and δn represents the displacement in normal direction. When particles move and rotate relative to each other, the contact point is modified. The normal component of the relative velocity at the contact point, denoted by ~unij , is given by: i  h ~unij = ~ui − ~u j · ~ni j ~ni j . (42) To calculate the normal contact force, the stiffness coefficient and the damping parameter in normal direction have to be determined. The normal stiffness coefficient can also be defined depending on the properties of the collision partners as: 4 kn = 3 r ri r j ri + r j  −1  1 − ν2 1 − ν2j  i  ,  + Ei Ej  (43) where νi , ν j and Ei , E j stand for the Poisson’s ratios and the Young’s moduli of the collision partners, respectively. The damping parameter in the normal direction, which represents the dissipation of energy during an inelastic collision, is determined by the following expression (Tsuji et al., 1992): q
ηn = αndam mi j kn δn 1/4 (44) with the empirical constant: αndam  2 ln en     − p 2 = π + ln2 en    2 for en , 0 (45) n for e = 0 where the symbol en denotes the coefficient of restitution in normal direction and is defined as the ratio of the normal relative velocities at the contact point before (index (0)) and after the contact: ~unij n e = ~uni j(0) = r h . h(0) (46) The modeling of the resulting tangential contact force is based on the spring-damper-slider system, where in this study the tangential contact force is assumed to be either a static or a sliding force (Link, 2006): t F~con,i j          t  d~ui j    = mi j =   dt           3/2 t~i j − ηt~uti j for F~ tstatic ≤ F~ tslide −kt δt | {z } F~ tstatic n ~ −µdyn F~con,i j ti j | {z } for F~ tstatic > F~ tslide , (47) F~ tslide where ~uti j stands for the tangential component of the relative motion at the contact:     ~uti j = ~ui − ~u j − ~unij + ri ω ~ j × ~ni j , ~ i + rj ω (48) 11 Powder Technology 264 (2014) 343–364 kt , ηt and δt represent the stiffness coefficient, the damping parameter and the displacement in the tangential direction, respectively, and µdyn denotes the dynamic friction coefficient. The damping parameter in the tangential direction ηt is defined by the expression: r  1/4 2 t t η = αdam m i j k t δt (49) 7 with the empirical constant:  2 ln et    − p   αtdam =  π2 + ln2 et    2 for et , 0 (50) for et = 0 where et is the restitution coefficient in the tangential direction. The tangential stiffness coefficient is calculated by using the normal stiffness coefficient and the material properties of the collision partners: (1 − νi )/Gi + (1 − ν j )/G j kt = . kn (1 − 0.5νi )/Gi + (1 − 0.5ν j )/G j (51) In order to calculate the time-dependent motion of the particle and fluid phases, different time steps are used in the DEM model, namely the fluid and the particle time steps. The size of the fluid time step can be arbitrarily selected and remains constant during the calculation. The particle time step is assumed to be constant as well in the DEM model. During each particle time step, the contact point is not modified and the particles do not move or rotate relative to each other (Alobaid, 2013). 2.2.3. Fluid-Particle interaction In this work, the employed program “DEMEST” is based on the DEM-code developed by Götz (2006) and extended by Alobaid et al. (2013). The program combines the classical computational fluid dynamics (CFD) to calculate the fluid phase with the discrete element method (DEM) to describe the solid phase. The balance equations of the single phase have to be extended in order to enable the fluid phase calculation in a gas-solid flow. For incompressible flows, the conservation equations of the gas phase can be written as: ∂ (αg ρg ) + ∇ · (αg ρg~ug ) = 0 ∂tg ∂ (αg ρg~ug ) + ∇ · (αg ρg~ug~ug ) = − αg ∇p + ∇ · (αg τg ) ∂tg (52) + αg ρg~g − β (~ug − ~u p ) Here ρg is the gas density, ~ug presents the velocity, p is the pressure, ~g stands for the gravitational acceleration, τg and β are the stress tensor for Newtonian fluids and the gas-particle momentum exchange coefficient, respectively. In the above equations, αg stands for the gas porosity (volume fraction) and is defined as: cv Np 1 X cv cv cv αg = 1 − α p = 1 − cv f V p,i . Vg i=1 p,i (53) where Vgcv is the volume of the fluid cell, V p,i stands for the volume of the particle i and N pcv represents the number of cv indicates the proportion of the particle volume, which exists in the control the particles in the fluid cell. The value f p,i volume. The momentum transfer is defined by the change of the velocities of the particles along their trajectories for each control volume based on all particles crossing this volume during a time interval between two consecutive fluid time steps. In the present work, a semi-implicit method is used to determine the momentum exchange according to Link (2006):  N cv  p X    ~ = β  ~ug − ~u p,i  , (54) R  i=1 12 Powder Technology 264 (2014) 343–364 where ~ug is the fluid velocity of the considered control volume, ~u p,i is the velocity of the particle i located in the control volume Vgcv . The symbol β denotes the inter-phase momentum transfer coefficient (resistance coefficient), which can be modeled using different approaches. However, the application of a certain model is only possible for a defined range of the porosity and the particle Reynolds number. Therefore, the selection of a model requires that both the porosity and the particle Reynolds number have to be known in the control volume. Various models to determine the resistance coefficient can be found in the literatures (e.g., Deen et al., 2007). In this work, the drag model by Hill et al. (2001) is applied, where the inter-phase momentum transfer coefficient for Re p > 40 is given by the following equation:    µg  (1 − αg )2 + B (1 − αg ) Re p  , β = 2 A αg dp (55) with the coefficients:   180, for αg < 0.6          A = 1  3 135   18 α3g 1 + √2 (1 − αg ) 2 + 64 (1 − αg ) ln(1 − αg ) + 16.14 (1 − αg )     , for αg ≥ 0.6   (1 − αg ) 1 + 0.681 (1 − α ) − 8.48 (1 − α )2 + 8.16 (1 − α )3 g g g B = 0.6057 (1 − αg )2 + 1.908 (1 − αg ) α2g + 0.209 α−3 g and the particle Reynolds number defined as: Re p = ρg d p ~u f − ~u p µg . (56) 13 Powder Technology 264 (2014) 343–364 3. Experiment and Numerical Set-up 3.1. Test rig 1.0 m The experiments have been carried out in a lab-scale cold-flow spouted fluidized bed model at two mass flow rates (two different fluidization conditions). Figure 2 shows the geometrical details of the experimental apparatus, which has a height of 1.0 m, a width of 0.15 m and a depth of 0.02 m. To allow visual observation of the particle motion inside the bed, the walls of the model are constructed of acrylic glass plates [Plexiglasr ]. m 0.007 m 0.17 m 0.02 0.001 m 0.15 m 0.008 m Figure 2: The spouted fluidized bed model used for the experiments. To obtain a homogeneous gas inflow and to suppress pressure drop fluctuations at the bottom of the bed, centrally located nine holes with 0.002 m inner diameter are used as gas distributors through which the fluidizing medium (air) is injected into the bed. The inlet flow rate is accurately controlled by mass flow controllers. To hinder the build-up of static electricity, the humidity of the injected air is maintained at 70%. The outlet located at the upper end of the apparatus is completely open to ensure an undisturbed air discharge. To prevent the particles from accidentally leaving the system, a coarse grid is mounted on the model outlet. The pressure drop over this grid is negligible. The solid phase consists of 36,500 spherical glass beads (particles) with 0.0025 m average diameter and a density of 2500 kg/m3 . The freeboard pressure is atmospheric, and the experiments were conducted at room temperature. The fluidizing air mass flow rates are chosen to be 0.005 kg/s and 0.006 kg/s, so that the glass particles do not leave the apparatus. Further parameters of both experiments are summarized in Table 1. In this experiment, the depth of the apparatus is assumed to be sufficiently small to achieve a pseudo-2D behavior which is necessary for recording the images using a high-resolution camera. The total measuring time lasted 500 ms. At both air mass flow rates, a series of snapshots of the bed behavior has been taken during the total measuring time for every 1 ms. 14 Powder Technology 264 (2014) 343–364 Gas Phase Density Dynamic viscosity Temperature Mass Flow Rates 1.225 kg/m3 1.79 × 10−5 kg/(m · s) 25 ◦ C 0.005 kg/s & 0.006 kg/s Solid Phase Geldart’s particle classification Total particle weight Density Average diameter Initial bed height D (Geldart, 1973) 0.75 kg 2500 kg/m3 0.0025 m 0.17 m Table 1: Parameters of the experimental set-up. Th high-speed camera used to capture the spatial distribution of particles in the present fluidized bed has the following properties: • Speed: 500 fps at full resolution • Maximum resolution: 1280 × 1024 pixel • Total recording time up to 2 seconds • Monochrome recording 3.2. Geometry and numerical grid The geometric model and the computational grid were created in ANSYSr WorkbenchTM 14.0 platform as depicted in Figure 3. Since the modeling of the inlet with nine holes would result in an unstructured mesh, unlike the experiment the inlet for the TFM simulations is modeled as a set of nine equivalent squares. To maintain the same mass flow rate of the injected gas and hence the same velocity, the area of a square is set equal to the area of a hole as shown on the left-side of the Figure 3. In this case, a block-structured grid with 80 × 216 × 22 cells was generated. A grid independence study was carried out, showing that a further grid refinement affects the solution very slightly (i.e., grid-independent solution). In the DEM simulations, the inlet is modeled as an equivalent nozzle (the area of the nozzle equals the area of the nine holes) due to restrictions on the used in-house code DEMEST as described in a previous work by Alobaid and Epple (2013). The reason for this assumption is that each hole of the inlet requires a separated sub-block to be generated. Since the inner diameter of each hole of the inlet is 0.002 m and the particle diameter is 0.0025 m, the generated sub-blocks in x and z directions will be smaller than the particle diameter. Due to the fact that the coupling between DEM and CFD becomes difficult when the CFD cells are smaller than the DEM particles, an additional grid, the so-called particle grid is implemented in the DEMEST code. This procedure allows the refinement of the fluid grid resolution beyond the particle scale (Alobaid et al., 2013). To minimize the computation time, the calculation of the DEM was parallelized, and therefore the computational domain is decomposed into eight similar blocks (1 2 3 4 & 6 7 8 9) and a central narrow block (5) for the inlet as displayed on the right-side in Figure 3. Note that the same number of cells is used for the DEM simulation to calculate the gas phase. 15 Powder Technology 264 (2014) 343–364 Outflow TFM Outflow DEM 1 2 3 4 56 7 8 9 0.007 m 0.008 m 0.007 m 0.008 m (experiment) (experiment) 0.02 m 0.007 m 0.008 m (simulation) (simulation) y z y x z Inflow x Domain Decomposition DEM Inflow Used Mesh TFM & DEM Figure 3: Numerical grid used for the simulations of both models. 3.3. Boundary and initial conditions Pressure measurements Figure 4 shows the boundary and initial conditions used for the simulations of both techniques (i.e., two-fluid model and discrete element method). The front, back, left- and right-side walls of the fluidized bed were modeled as Outflow located at the bottom of the bed was impermeable which means no-slipOutflow rigid walls for both phases. The distributor 1 2 3 4 5 6 At 7 8 9the inlet, an uniform conceived as a prescribed gas phase and an impermeable no-slip rigid wall for the solid phase. distribution of the velocity components of the gas phase is assumed. At the initial time, the velocity of the injected gas through the inlet (the nine holes for the TFM or the nozzle for the DEM simulations) was increased instantaneously from minimum fluidization velocity to the required spout velocity. Only two air mass flow rates, 0.005 kg/s and 0.006 kg/s , are considered in this study. The pressurized air enters the fluidized bed at these mass flow rates, and the solids mass flow rate is set equal to zero. At the top of the bed, the outlet is a continuous outflow for the fluid phase. The pressure is P4 set at an ambient atmosphere. For the two-fluid model simulations, the wall boundary condition by Johnson and Jackson (1987) is applied for the tangential velocity and the granular temperature ofDEM the solid phase at the wall. As depicted in TFM P Figure 4, a freeboard of30.83 m was provided to allow for the bed expansion. In accordance withP the experiments, the minimum fluidization conditions were prescribed as the initial conditions 2 for the numerical calculations of both techniques which implies a force balance between the buoyant weight of the sus1 the forces exerted by the fluidizing gas. As a consequence of this hydrodynamic equilibrium, pended solid particles Pand no net particle movement occurs. As displayed in Figure 4, the initial bed height was identical to the experimental minimum fluidization bed height (0.17 m) as explained in Section 3.1. For the two-fluid model simulations, the particles y y are lying in the spouted bed with a maximum packing limit (αmax = z0.63). The initial velocity field and pressure of the z Domain Decomposition x Fluidized Bed Model Used Mesh gas and solids phases are set to zero. The inlet massTFM flow rate of the gasx phase to 0.005 kg/s for the first Inflow is set equalDEM Inflow EXP. & DEM case and 0.006 kg/s for the second one, while the mass flow rate of the solid phase for both cases is set to zero. The turbulence is specified at the inlet by the turbulent intensity (10%) and the hydraulic diameter. 16 . mg = 0.006 [kg/s] Powder Technology 264 (2014) 343–364 pgauge=0 [Pa] pgauge=0 [Pa] inlet wall outlet pgauge=0 [Pa] αg= 1 αs= 0 pg= 0 ug= 0 Freeboard H=1.0 m Freeboard αg= 1 αp= 0 pg= 0 ug= 0 H0=0.17 m αg= 0.37 αs= 0.63 ug= 0 . mg,1 = 0.005 [kg/s] . mg,2 = 0.006 [kg/s] . mg = 0.005 & 0.006 [kg/s] . ms = 0 [kg/s] . mg = 0.005 & 0.006 [kg/s] . ms = 0 [kg/s] Experiments ααgg==0.37 0.37 ααps= 0.63 =0.63 ug= 0 TFM DEM Figure 4: Boundary and initial conditions for the simulations of both techniques. 3.4. Simulation settings The simulations were carried out in diverse settings such as drag models, particle wall boundary conditions and solver settings (the discretization scheme, time step, etc.). The commercial CFD package ANSYSr FLUENTTM 14.0 was employed for the Euler-Euler (two-fluid model) simulations, and the in-house code DEMEST (Alobaid et al., 2013) was used for the Euler-Lagrange (discrete element method) simulations. The corresponding numerical settings for the two-fluid model and the discrete element method simulations are listed in Table B.2 and Table C.3, respectively. In both models, the properties of the gas phase are the same, but the parameters of the solid phase are dependent on the model. Note that the input parameters of the two-fluid model listed in Table B.2 are based on the comparative study done in Section 4.1. Furthermore, the settings used in the DEM simulation, listed in Table C.3, were selected on the basis of the studies of Kharaz et al. (1999, 2001) and Di Renzo and Di Maio (2004), so that these settings result in the most appropriate results. Here, it is worth noting that the stiffness coefficient is a crucial parameter for obtaining correct physical results. A stiffness coefficient with a high order of magnitude (106 − 108 N/m) requires very small time steps, which can hardly be achieved with currently available computing performance. However, in these CFD/DEM simulations, the stiffness coefficient is set equal 4.1 × 105 N/m since it shows a very good compromise between the acceptable computing time and the achieved accuracy. This value was already derived through a parameter study by Alobaid et al. (2014). 4. Results and Discussion In this work, two different modeling strategies for the simulation of the particulate flow in a spouted fluidized bed are adopted, the Euler-Euler approach (two-fluid model) and the Euler-Lagrange approach (discrete element method). In this section, the numerical results of both techniques are compared against each other and with the experimental data of particulate flows in a cold-flow lab-scale gas-solid spouted fluidized bed reactor at two different operating conditions. The most critical comparison between experiments and the simulation results using both techniques is given by the analysis of the bubble size and the bed expansion. To evaluate the process of bubble formation quantitatively, the 17 Powder Technology 264 (2014) 343–364 circular equivalent bubble diameter, denoted here by deq , is introduced: r 4 Ab , deq = π (57) where Ab stands for the measured bubble area calculated from the experimental snapshots and the computational contours of the solids volume fraction. 4.1. Influence of parameters of the two-fluid model To obtain the best simulation results, the influence of various simulation parameters such as drag models, solidphase wall boundary conditions in terms of the specularity coefficient, the particle-particle restitution coefficient and the granular temperature approach on the global bed hydrodynamics has been investigated. Here, the parameters and models were selected based on the results and recommendations of previous works mentioned in Section 1.1 and 1.2. 4.1.1. Gas-solid drag model In gas-solid flows, the momentum exchange between both phases characterized by a drag force is caused by the relative motion between the phases. This force whose formulation is based on the force acting on a single particle in a fluid flow has been well investigated and empirically correlated to a wide range of particle Reynolds numbers. Since there are several correlations reported in the literature for calculating the gas-solid momentum exchange coefficient, the appropriate choice of a drag law is required. In the present work, three drag models were tested, namely Syamlal and O’Brien (1989), Gidaspow et al. (1992) and Wen and Yu (1966). For these tests, the particle-particle restitution and specularity coefficients are set to 0.6 and 0.2, respectively, and the algebraic granular temperature approach is applied. The bed height and the equivalent bubble diameter obtained at both mass flow rates as a function of the time are shown in Figure 5. 50 Exp. Gidaspow et al. Syamlal-O'Brien Wen-Yu Exp. Gidaspow et al. Syamlal-O'Brien Wen-Yu 45 Bed height [cm] Bed height [cm] 30 25 20 40 35 30 25 20 15 15 0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500 Time [ms] Time [ms] 25 Exp. Gidaspow et al. Syamlal-O'Brien Wen-Yu Equivalent diameter [cm] Equivalent diameter [cm] 15 10 5 Exp. Gidaspow et al. Syamlal-O'Brien Wen-Yu 20 15 10 5 0 0 0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500 Time [ms] Time [ms] Figure 5: Influence of drag models on the bed behavior at a gas mass flow rates of 0.005 kg/s (left) and 0.006 kg/s (right). According to the above diagrams, it can be deduced that the most appropriate drag model for the present set-up is that of Syamlal and O’Brien (1989). This drag model reproduces the flow pattern (bed expansion and bubble formation) better than the other tested models for both mass flow rates. In the results of the Wen and Yu (1966) drag model, it is 18 Powder Technology 264 (2014) 343–364 clear that the bed height starts to decrease at 300 ms, however, the equivalent diameter does not change the increasing trend until 340 ms, which is a discrepancy itself. When compared to the experimental bed behavior, a small deviation in the predicted bed expansion and the bubble formation using Syamlal and O’Brien (1989) drag model during the final stage at a gas mass flow rate 0.005 kg/s can be noticed, and a relatively large deviation of the predicted bed behavior at a gas mass flow rate 0.006 kg/s in the range of time between 75 ms and 500 ms. This discrepancy shall be reduced by further investigating the effect of other simulation parameters on the flow behavior which is the topic of the next subsections. 4.1.2. Restitution coefficient As described before, the formulation of the drag force is based on the force acting on a single particle in a fluid flow. However, when a single particle moves in a dispersed two-phase mixture, the drag is affected by the presence of other particles, and therefore the influence of the particle-particle interaction must be taken into account. The elasticity of the particle-particle collision can be modeled by the restitution coefficient. In the previous step, the drag model by Syamlal and O’Brien (1989) was chosen as the best drag model among the tested drag models. In this step, it is adapted and the influence of four different values (0.60, 0.70, 0.80 and 0.90) of the particle-particle restitution coefficients on the hydrodynamics of the investigated spouted bed is studied. In addition, the specularity coefficient of 0.2 is chosen, and the algebraic granular temperature approach is used. The outcome of this parameter study is depicted in Figure 6, which again shows the bed height and equivalent bubble diameter as a function of the time. 50 Exp. RC = 0.6 RC = 0.7 RC = 0.8 RC = 0.9 25 Exp. RC = 0.6 RC = 0.7 RC = 0.8 RC = 0.9 45 Bed height [cm] Bed height [cm] 30 20 40 35 30 25 20 15 15 0 50 0 100 150 200 250 300 350 400 450 500 50 100 150 200 250 300 350 400 450 500 Time [ms] Time [ms] 30 Exp. RC = 0.6 RC = 0.7 RC = 0.8 RC = 0.9 10 Equivalent diameter [cm] Equivalent diameter [cm] 15 5 Exp. RC = 0.6 RC = 0.7 RC = 0.8 RC = 0.9 25 20 15 10 5 0 0 0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500 Time [ms] Time [ms] Figure 6: Influence of the particle-particle restitution coefficient (RC) on the bed behavior at a gas mass flow rates of 0.005 kg/s (left) and 0.006 kg/s (right). Looking at Figure 6, it is clear that the increase of the particle-particle restitution coefficient results in different flow structures (lower bed heights and smaller bubble sizes) at the same instant of time. The overall assessment of the effect of the particle-particle restitution coefficient on the computed bed behavior shows that the calculation with a restitution coefficient of 0.60 yields a good agreement with the experimental data. However, there is still a large deviation between the predicted bed behavior and the experimentally observed one at a gas mass flow rate of 0.006 kg/s. It should also be noticed that, after the bubble started to burst at approximately 350 ms for the gas mass flow rate of 0.005 kg/s, the computed bed height and bubble size become smaller than the experimental ones. This implies that the predicted solids layer above the bubble is thinner than the experimental one, and hence the bubble burst starts a bit earlier than observed 19 Powder Technology 264 (2014) 343–364 in the experiment. It is worth noting that the use of a restitution coefficient smaller than 0.6 resulted in overestimating both the bed height and the bubble size at the gas mass flow rate of 0.005 kg/s. Furthermore, the baseline value of the restitution coefficient 0.6 is kept for the simulations when the gas mass flow rate is increased to 0.006 kg/s, since our goal is to use the same model settings. 4.1.3. Granular temperature approach In the previous steps, the algebraic granular temperature approach has been adapted while testing the influence of various parameters of the two-fluid model. In this step, the best drag model and the restitution coefficient are taken and the significance of the partial differential approach of the granular temperature is revealed. As depicted in Figure 7 showing again the same quantities as before, the partial differential approach overpredicts the bed expansion and the bubble sizes at the gas mass flow rate of 0.005 kg/s in the time interval between 0 and 250 ms and underpredicts them during the remaining time of the process of bubble formation. Although both approaches underpredict the bed behavior at the gas mass flow rate of 0.006 kg/s, the algebraic model gives much better results than the partial differential one. It is noteworthy that the algebraic approach for the granular temperature requires a lower computational effort than the partial differential model, since the convection and the diffusion terms in the transport equation of the granular temperature are not taken into account. Therefore, the algebraic approach of the granular temperature is used by default for the next steps. In the previous steps, the chosen simulation parameters could somehow predict the flow behavior. However, to further improve the computational results, we shall look in the next step at the influence of the particle-wall boundary conditions in terms of the specularity coefficient on the hydrodynamic behavior of the investigated gas-solid flows. 30 50 Exp. PDE Algebraic 45 Bed height [cm] Bed height [cm] Exp. PDE Algebraic 25 20 40 35 30 25 20 15 15 0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500 Time [ms] Time [ms] 25 Exp. PDE Algebraic Equivalent diameter [cm] Equivalent diameter [cm] 15 10 5 0 Exp. PDE Algebraic 20 15 10 5 0 0 50 100 150 200 250 300 350 400 450 500 Time [ms] 0 50 100 150 200 250 300 350 400 450 500 Time [ms] Figure 7: Influence of the granular temperature approach on the bed behavior at a gas mass flow rates of 0.005 kg/s (left) and 0.006 kg/s (right). 4.1.4. Particle-wall boundary conditions In this step, the influence of the particle-wall boundary condition in terms of specularity coefficient ϕ on the hydrodynamics of the particulate flow is examined. The specularity coefficient is an empirical parameter used to quantify the nature of the particle-wall collisions, and it characterizes the friction between particles and the wall. The values of the specularity coefficient have to be in the interval [0,1]. A value of zero implies a specular collision (free-slip conditions) during which there is no friction between the particles and the wall, and hence the shear stress is zero. A value of 20 Powder Technology 264 (2014) 343–364 unity corresponds to a diffuse collision (no-slip conditions) during which the particles stick to the wall, and thus the shear stress is the maximum. To evaluate the impact of the specularity coefficient on the flow behavior in the present spouted bed and to obtain the most suitable value, five different specularity coefficients (ϕ = 0, 0.2, 0.5, 0.7 and 1.0) are investigated. As shown in Figure 8, it is clear that the specularity coefficient has a pronounced effect on both bubble size and bed height. It is also obvious that the bed height is the highest for ϕ = 0, since there is no friction between the particles and the wall. This free-slip of particles on the wall increases their velocity and consequently results in a higher bed height. On the contrary, the no-slip boundary condition (i.e., ϕ = 1) results in a low bed height, because a strong friction exists between the particles and the wall, and hence it increases the resistance of the downward flow of particles. Based on the various factors considered and on the basis of the information available, it can be concluded that the value of ϕ = 0.5 provides the best agreement with the experimental bed behavior at a gas mass flow rate of 0.005 kg/s. 50 Exp. SC= 0.0 SC= 0.2 SC= 0.5 SC= 0.7 SC= 1.0 25 Exp. SC = 0.0 SC = 0.2 SC = 0.5 SC = 0.7 SC = 1.0 45 Bed height [cm] Bed height [cm] 30 20 40 35 30 25 20 15 15 0 50 0 100 150 200 250 300 350 400 450 500 50 100 150 200 250 300 350 400 450 500 Time [ms] Time [ms] 30 Exp. SC= 0.0 SC= 0.2 SC= 0.5 SC= 0.7 SC= 1.0 15 Equivalent diameter [cm] Equivalent diameter [cm] 20 10 5 0 Exp. SC = 0.0 SC = 0.2 SC = 0.5 SC = 0.7 SC = 1.0 25 20 15 10 5 0 0 50 100 150 200 250 300 350 400 450 500 Time [ms] 0 50 100 150 200 250 300 350 400 450 500 Time [ms] Figure 8: Influence of the particle-wall boundary conditions in terms of the specularity coefficient (SC) on the bed behavior at a gas mass flow rates of 0.005 kg/s (left) and 0.006 kg/s (right). 4.1.5. Summary of key findings The previous parameter study examined several settings of the two-fluid model including the three drag models, solid-phase wall boundary conditions in terms of the specularity coefficient, the particle-particle restitution coefficient and the granular temperature approach. For the present spouted fluidized bed, the obtained results at the gas mass flow rate of 0.005 kg/s provide evidence of the following findings: • The drag model has a significant effect on the predicted flow behavior. The drag model of Syamlal and O’Brien (1989) is the most suitable one to investigate the hydrodynamics of the present gas-solid fluidized bed. Note that, in this step, the baseline value of the specularity coefficient (ϕ = 0.2) was selected based on Lan et al. (2012). Using these settings (i.e., the drag model and the value of the specularity coefficient), different values of the particle-particle restitution coefficient were tested (e ss = 0.9, 0.8, 0.7 and 0.6). A value of 0.6 for the particleparticle restitution coefficient gives the best agreement between the predictions of the two-fluid model and the experiment. • The algebraic granular temperature approach is more appropriate than the partial differential model, since it reproduces the flow structure reasonably and is more efficient than the partial differential approach. 21 Powder Technology 264 (2014) 343–364 • The use of a suitable wall boundary condition is critical for improving the predictions of the two-fluid model. The impact of the specularity coefficient on the bed behavior has been evaluated using different values. In the range of values considered for the specularity coefficient (ϕ = 0, 0.2, 0.5, 0.7 and 1.0), a value of 0.5 reproduces the flow patterns at best. 4.2. Comparison of Euler-Euler and Euler-Lagrange models Once the steps to “optimize” the results have been established in the preceding section, the attention is now turned to the evaluation of the overall bed behavior at both gas mass flow rates. Here, the numerical simulation results obtained by the two-fluid model and the discrete element method are compared against the experiment. A critical comparison of the bed expansion is first presented. Then, the process of the bubble formation is evaluated in terms of the equivalent circular bubble diameter. 4.2.1. Gas mass flow rate of 0.005 kg/s To validate the numerical results obtained in the framework of the two-fluid model and the discrete element method, a critical comparison of the predicted results against the experiment at the gas mass flow rate of 0.005 kg/s is displayed in Figure 9 and 10. Here, the experimental snapshots, the distribution of the vertical component of the particle velocity obtained by the discrete element method and the numerical contours of the solids volume fraction using the two-fluid model are plotted, so that the most important features of the bubble growth and the bed height are shown. At the outset of the experiment, the initial bed height is identical to the experimental minimum fluidization bed height. Once a pressurized air enters the fluidized bed at the mass flow rate of 0.005 kg/s, the particles become suspended and a small bubble is formed at the inlet. As the bubble is growing with the time, the layer of particles above the bubble becomes thinner and moves faster. After approximately 375 ms, the bubble starts to burst, then the particles fall downwards to the spout region due to the gravitational force. In the time period between 0 and 300 ms, the numerical results of the two-fluid model agree almost perfectly with the experimental data. Beyond 300 ms, the two-fluid hydrodynamic model reproduces neither the experimental bed expansion nor the bubble sizes in a close agreement with the experiment, since it tends to predict somewhat smaller bubbles and lower bed heights. It is also worth to notice that this model still shows a high correlation concerning the experimental jet region during the final stage of the process of the bubble formation. In that period of time, the top layer of the upwards moving particles loses more and more of its mass, and the particles fall down due to the gravitational force. Hence, the predicted bubble starts to burst earlier than in the experiment. On the contrary, the discrete element method could not predict the bubble shape with similar accuracy as the two-fluid model during the whole time, except for the jet region. However, for the present test case, the predictions of the DEM are in reasonable agreement with the experiment. 22 Powder Technology 264 (2014) 343–364 [m] 0.5 0.4 0.3 0.2 0.1 0.0 0.5 0.4 0.3 0.2 0.1 0.0 -1 0 +1 vp [m/s] 0.5 0.4 0.3 0.2 0.1 0.0 0 [ms] 25 [ms] 50 [ms] 75 [ms] 0 0.1 100 [ms] 0.2 125 [ms] 0.3 0.4 150 [ms] 0.5 175 [ms] 0.6 200 [ms] 225 [ms] αs [-] Figure 9: Snapshots of the solid phase at a gas mass flow rate of 0.005 kg/s in the time range between 0 and 225 ms: (upper row) experiment, (middle row) vertical component of the particle velocity predicted by the DEM and (lower row) volume fraction predicted by the TFM. 23 Powder Technology 264 (2014) 343–364 [m] 0.5 0.4 0.3 0.2 0.1 0.0 0.5 0.4 0.3 0.2 0.1 0.0 -1 0 +1 vp [m/s] 0.5 0.4 0.3 0.2 0.1 0.0 250 [ms] 275 [ms] 300 [ms] 325 [ms] 0 0.1 350 [ms] 0.2 375 [ms] 0.3 0.4 400 [ms] 0.5 425 [ms] 0.6 450 [ms] 475 [ms] αs [-] Figure 10: Snapshots of the solid phase at a gas mass flow rate of 0.005 kg/s in the time range between 250 and 475 ms: (upper row) experiment, (middle row) vertical component of the particle velocity predicted by the DEM and (lower row) volume fraction predicted by the TFM. Figure 11 depicts the experimentally observed and the computed bed behavior as a function of time. It shows a good agreement between the bed expansion and bubble formation predicted by the two-fluid model and the experimental data. In the time interval between 0 and 340 ms, the two-fluid model predicts the experimental bed heights in a nearly perfect agreement. Furthermore, the diagrams indicate that the two-fluid model tends to predict smaller expansion and bubble sizes than the experimentally observed ones during the final stage of the bubble formation, especially as the bubble starts to burst. On the other hand, the discrete element method could predict neither the bed expansion nor the bubble size during the whole process with the same level of agreement as the two-fluid model. 24 Powder Technology 264 (2014) 343–364 15 30 Equivalent diameter [cm] Bed height [cm] Exp. EL/DEM EE/TFM 25 20 Exp. EL/DEM EE/TFM 10 5 0 15 0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500 Time [ms] Time [ms] Figure 11: Experimental and computed bed expansion (left) and bubble growth (right) using both simulation techniques as a function of time at a gas mass flow rate 0.005 kg/s. 4 Vertical particle velocity [m/s] Vertical particle velocity [m/s] To make a further judgment about the success of the simulations, the vertical particle velocity profiles on the midplane (i.e., y=0.01 m) at different locations along the bed height and instants of time are plotted in Figure 12. The selected positions for the comparison ensure that the three regions of the spouted bed (i.e., spout, annulus and fountain) are covered. When looking at the line located at 0.12 m, it is clear that, in the spout region, the particles rise at high velocities, while they move slowly downwards in the annulus region between the spout and the walls. At the instant of time 350 ms, the negative values of the vertical particle velocities along the line located at 0.2 m refer to a strong downward movement of the particles in the annulus region. At the same time, the particles in the fountain region move slowly upwards along the line located in the fountain region. This behavior changes at 450 ms after the bubble has already burst. At this time, the particles in the annulus region are circulating, while the particles in the fountain region are separated towards the apparatus wall, and afterwards they move back to the jet zone due to the gravity force. This behavior can be seen in Figure 10 during the final stage. It is worth mentioning that the vertical velocity of the particles predicted by the two-fluid model are not zero, since the specularity coefficient is selected less than unity. Thus, the particles do not stick to the wall as described in Section 4.1.4. Further inspections of Figure 12 shows a very good agreement between the predictions of both simulation techniques, except on the wall, where the DEM tends to predict somewhat larger values of the vertical component of the particle velocity. Note that, in this comparison, the experimental results are not plotted since they are not available due to the lack of the required measuring techniques. DEM: 50[ms] TFM: 50[ms] DEM: 250[ms] TFM: 250[ms] 3 2 1 0 -1 -2 0 2.5 5 7.5 10 12.5 4 2 1 0 -1 -2 15 x [cm] DEM: 350[ms] TFM: 350[ms] DEM: 450[ms] TFM: 450[ms] 3 0 2.5 5 7.5 10 12.5 15 x [cm] Figure 12: Vertical component of the particle velocity for a gas mass flow rate of 0.005 kg/s at different locations and instants of time. The measuring levels are located at a height of 0.12 m (left) and 0.20 m (right). According to the above discussed evaluation procedures including the qualitative comparison of the bed behavior shown in Figure 9 and 10 and the combination of the results in Figure 11 and 12, it can be concluded that the two-fluid model can predict both the bubble size and the bed expansion of the investigated spouted bed better than the EulerLagrange model at gas mass flow rate of 0.005 kg/s during the whole process of the bubble formation. Note that the assumption made at the inlet (a nozzle instead of nine holes) represents one of the most probable reasons for deviations 25 Powder Technology 264 (2014) 343–364 of the DEM model. 4.2.2. Gas mass flow rate of 0.006 kg/s As mentioned before, the successfully validated TFM settings and the used DEM simulation parameters for the particulate flow in the investigated fluidized bed model at a gas mass flow rate of 0.005 kg/s are maintained, while only the gas flow rate is increased to 0.006 kg/s. In this section, the numerical simulation results obtained by adopting the two-fluid model and the discrete element method are compared against the experimental data at the gas mass flow rate of 0.006 kg/s. In Figure 13 and 14, the computed contours of the spatial distribution of the vertical velocity component of the particles based on the discrete element method are plotted, while the contours of the solids volume fraction are given in the framework of the two-fluid model. Here, the same parameters of both models as chosen for the gas mass flow rate of 0.005 kg/s are used. It is noticed that the typical flow patterns of a spouted bed including three regions, the spout in the center, the fountain above the initial bed surface and the annulus between the spout and the walls, are clearly observed using both simulation techniques. 26 Powder Technology 264 (2014) 343–364 [m] 0.5 0.4 0.3 0.2 0.1 0.0 0.5 0.4 0.3 0.2 0.1 0.0 -1 0 +1 vp [m/s] 0.5 0.4 0.3 0.2 0.1 0.0 0 [ms] 25 [ms] 50 [ms] 75 [ms] 0 0.1 100 [ms] 0.2 125 [ms] 0.3 0.4 150 [ms] 0.5 175 [ms] 0.6 200 [ms] 225 [ms] αs [-] Figure 13: Snapshots of the solid phase at a gas mass flow rate of 0.006 kg/s in the time range between 0 and 225 ms: (upper row) experiment, (middle row) vertical component of the particle velocity predicted by the DEM and (lower row) volume fraction predicted by the TFM. 27 Powder Technology 264 (2014) 343–364 [m] 0.5 0.4 0.3 0.2 0.1 0.0 0.5 0.4 0.3 0.2 0.1 0.0 -1 0 +1 vp [m/s] 0.5 0.4 0.3 0.2 0.1 0.0 250 [ms] 275 [ms] 300 [ms] 325 [ms] 0 0.1 350 [ms] 0.2 375 [ms] 0.3 0.4 400 [ms] 0.5 425 [ms] 0.6 450 [ms] 475 [ms] αs [-] Figure 14: Snapshots of the solid phase at a gas mass flow rate of 0.006 kg/s in the time range between 250 and 475 ms: (upper row) experiment, (middle row) vertical component of the particle velocity predicted by the DEM and (lower row) volume fraction predicted by the TFM. At the beginning, the initial bed height is identical to the experimental minimum fluidization bed height. Once the pressurized air enters the fluidized bed at the mass flow rate 0.006 kg/s, the particles become suspended and a small bubble is formed. As the bubble is growing with the time, the layer of particles above the bubble becomes thinner and moves faster as visible in the experimental snapshots in Figure 13 and 14. After approximately 350 ms, the bubble starts to burst in the experiment, and then the particles fall downwards to the spout region due to the gravitational force. In the time interval between 0 and 75 ms, the two-fluid model reproduces the ellipsoidal shape of the experimental bubble very reasonably, while a circular one is predicted by the discrete element model. However, it is clearly obvious that, during this time, the DEM results appear to underestimate the jetting of particles into the bubble, while the TFM results overestimate the distance traveled by this jet. Beyond this period of time, the two-fluid model tends to predict somewhat smaller bubbles and lower bed heights, and hence it predicts neither the experimental bed expansion nor the bubble sizes. 28 Powder Technology 264 (2014) 343–364 In the time interval between 0 and approximately 225 ms, the discrete element model underpredicts the experimental bed behavior, but the bubble shape is similar to the experimental one. It is also noticed that the TFM still shows high correlation towards the experimental jet, while a small deviation occurs in the DEM predictions and becomes larger during the remaining time of the process of the bubble formation. In the time interval between 250 and 350 ms, the top layer of the upwards moving particles loses more and more of its mass, and the particles fall down along the side walls due to the gravitational force. In this period of time, the DEM predicts the bed expansion in close agreement to the experiment while the deviation in the two-fluid model results becomes larger than before. During the final stage of the bubble formation, the discrete element method predicts the experimental bed expansion very reasonably, while the maximum deviation in the two-fluid model predictions is reached. In this period of time, the two-fluid model can only reproduce the spout region, while unrealistic flow behavior is noticed in the fountain region. Figure 15 displays the experimental and the computed dynamic bed height and bubble size using both techniques as a function of time. A good agreement between the predictions of the discrete element method and the experimental data is shown. A large deviation between the two-fluid model results and the experiment is observed from 75 ms till the end of the observed time of bubble formation. 25 50 Bed height [cm] 45 Equivalent diameter [cm] Exp. EL/DEM EE/TFM 40 35 30 25 20 Exp. EL/DEM EE/TFM 20 15 10 5 0 15 0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500 Time [ms] Time [ms] Figure 15: Experimental and computed bed expansion (left) and bubble growth (right) using both simulation techniques as a function of time at a gas mass flow rate 0.006 kg/s. Analog to the comparison done in the preceding section, the vertical particle velocity profiles on the midplane at different locations along the bed height and instants of time are plotted in Figure 16. When looking at the line located at 0.15 m, it is clear that, in the spout region, the particles rise at high velocities, while they move slowly downwards in the annulus region between the spout and the walls. At the instant of time 50 ms, the predictions of both models agree well, while a large deviation between them can be observed at 250 ms which confirms the qualitative predictions shown in Figure 13 and 14 and the macroscopic outputs plotted in Figure 15. At the instants of time 350 ms and 450 ms, the negative values of the vertical particle velocities along the line located at 0.25 m refer to a strong downward movement of the particles (circulation) in the annulus region. Note that the two-fluid model predicts zero velocity in the region near to the wall since no particles are there which can obviously be observed in the contours of the solids volume fraction demonstrated in Figure 14. At the instant of time 450 ms, the predicted velocities are smaller than those at 350 ms, since the bubble has already burst. Thus, the particles in the annulus region are circulating, while the particles in the fountain region are separated towards the apparatus wall, and afterwards they move back to the jet zone due to the gravity force. This behavior can be seen in Figure 14 during the final stage. Overall, the discrete element model can capture the experimental flow behavior better than the two-fluid model at the gas mass flow rate of 0.006 kg/s. The possible reasons for the deviation are explained in Section 4.3. 29 4 Vertical particle velocity [m/s] Vertical particle velocity [m/s] Powder Technology 264 (2014) 343–364 DEM: 50[ms] TFM: 50[ms] DEM: 250[ms] TFM: 50[ms] 3 2 1 0 -1 -2 0 2.5 5 7.5 10 12.5 4 2 1 0 -1 -2 15 x [cm] DEM: 350[ms] TFM: 350[ms] DEM: 450[ms] TFM: 450[ms] 3 0 2.5 5 7.5 10 12.5 15 x [cm] Figure 16: Vertical component of the particle velocity for a gas mass flow rate of 0.006 kg/s at different locations and instants of time. The measuring levels are located at a height of 0.15 m (left) and 0.25 m (right). 4.3. Discussion The overall assessment of the numerical results obtained by the two-fluid model shows that this model can reproduce appropriately the experimental flow pattern at the gas mass flow rate of 0.005 kg/s, but it does not succeed in reproducing them when the gas mass flow rate is increased to 0.006 kg/s, especially during the final stage of the bubble formation. This can probably be attributed to the fact that TFM does not recognize the discrete nature of the solid phase, which effects the calculation of the particle-particle interactions and the drag force (Gera et al., 1998). Furthermore, the assumptions made in the selected drag model are insufficient when increasing the gas mass flow rate to 0.006 kg/s as well as the consequence of neglecting the particle rotation. At the gas mass flow rate of 0.006 kg/s, the numerical results based on the discrete element method agree well with the experiment, except in the spout region during the period of time from 250 ms to 500 ms. These differences can be attributed to several reasons: The most important one is related to the inaccurate solid-phase parameters used in this calculation, since the simulation of the solids phase in the framework of discrete particle models requires a lot of material and geometrical properties of the particle and wall pairing, which are not completely available. Another important reason is related to the gas-solid interactions in the spout region, where the maximum velocity difference exists between the phases. The empirical correlations applied to calculate the drag coefficient produce an extra error concerning the determination of the momentum exchange, in particular when the velocity differences between the phases are high (Bokkers et al., 2004). Another reason for the discrepancy of the numerical model is related to the absence of turbulence modeling in the DEM simulations. These errors in the calculation of the fluid phase are forwarded to the particle phase by means of data transfer between both phases, and the particle calculation with DEM is associated itself with numerical errors such as the penetration depth between the collision partners, especially if the stiffness coefficient is underestimated. Overall, using the same simulation parameters of each model for both operation conditions, the two-fluid model predicts the bed behavior at a gas mass flow rate of 0.005 kg/s more precisely than the discrete element model, while DEM could capture the complex flow structures much better than the two-fluid model at a gas mass flow rate of 0.006 kg/s. However, the DEM requires much longer computation time, because each individual particle must be tracked. Hence, a set of ordinary differential equations has to be solved for each particle which is one of the main drawbacks of the Euler-Lagrange model. Although the stiffness coefficient was reduced for the present simulations, the computation time for DEM was approximately two times longer than the time required for the two-fluid model. 5. Conclusions To gain insight into the hydrodynamics of a cold-flow quasi-2D spouted fluidized bed, a comparative study of two modeling approaches based on a global comparison with experimental observations for bubble size in terms of the representative bubble diameter and the bed height has been carried out. The main goal is to explore the effects of the flow rate of the injected gas on the ability of the used simulation techniques, the Two-Fluid Model and the Discrete Element Method, to predict the flow behavior in a lab-scale spouted fluidized bed for two different operating conditions while maintaining the same settings. To obtain the best simulation settings in the framework of TFM, the influence 30 Powder Technology 264 (2014) 343–364 of various simulation parameters such as gas-particle drag models, the particle wall boundary conditions in terms of the specularity coefficient, the particle-particle restitution coefficient and the granular temperature approach on the flow structure has been investigated. The model settings are validated and used for the simulation at a higher gas mass flow rate with the same computational set-up. Euler-Euler Model The overall assessment of the volume fraction contours of the solid phase predicted in the framework of the two-fluid model shows the following: • The drag model has a significant effect on the simulation results, hence on the flow behavior. The drag model of Syamlal and O’Brien (1989) is suitable to investigate the hydrodynamics of the present gas-solid fluidized bed. However, the agreement between numerical simulation and experiment is not perfect, especially at a gas mass flow rate of 0.006 kg/s, which can be attributed to the fact that the choice of an appropriate drag model alone might probably not improve the simulation results, and other contributions such as a frictional viscosity model, the particleparticle interactions and the particle rotation are not appropriately modeled or neglected in the framework of the two-fluid model. Note that the energy dissipation due to particle rotation is partially captured by the restitution coefficient (see Eq. (17)). • In the range of values considered for the particle-particle restitution coefficient (e ss = 0.6 to 0.9), a value of 0.6 gives the best predictions at a gas mass flow rate of 0.005 kg/s. Therefore, the same value is used for the simulations at both operating conditions. • The particle wall boundary conditions in terms of the specularity coefficient have a significant effect on the flow behavior of spouted beds characterized by the bed height and the bubble size. Thus, the use of a suitable boundary condition is critical for improving the predictions. The impact of the specularity coefficient on the bed behavior has been evaluated by testing different values. A value of ϕ = 0.5 reproduces the best predicted flow patterns. • Surprisingly, the algebraic granular temperature approach is more appropriate than the partial differential one at both mass flow rates, since it could reproduce the flow structure reasonably and is more efficient than the partial differential approach. The two-fluid model, despite its ability to capture the flow features obtained in the experiment at a mass flow rate of 0.005 kg/s, does not succeed in reproducing the experiment when increasing the mass flow rate to 0.006 kg/s, especially during the final stage of bubble formation. Therefore, to get more realistic results, additional investigations are required to overcome this shortcoming. Euler-Lagrange Model According to the simulation results based on the discrete element method, a good agreement between the numerical predictions and the experiments of the present lab-scale gas-solid fluidized bed at both mass flow rates can be obtained. Although the discrete element model could reproduce the flow pattern more appropriately than the two-fluid model at the mass flow rate of 0.006 kg/s, it requires a much longer computation time, since a set of ordinary differential equations is solved for each particle. However, DEM and TFM are methods, which do not target at the same process level. TFM is a good approach for macroscale simulations, whereas DEM is a tool for microscopic investigations, since it is restricted in both process time and scale to laboratory facilities. Furthermore, the lack of turbulence modeling in the CFD/DEM model is a serious shortcoming of this model. On the basis of the above considered factors, it can be concluded that, for the present spouted fluidized bed, the two-fluid model is able to simulate the complex flow in the present spouted fluidized bed at a gas mass flow rate of 0.005 kg/s better than the discrete element method, while the latter is more successful when increasing the gas mass flow rate to 0.006 kg/s. However, the predicted characteristics of the bubble formation and the bed expansion using DEM at a gas mass flow rate of 0.005 kg/s are in a good agreement with the experimental observation. When the errors produced by the measurement are considered, we think that the deviation in the bed expansion and the bubble formation is acceptable. 31 Powder Technology 264 (2014) 343–364 Appendix A. Force and Momentum Models for DEM The considered fluid forces acting on a particle are given by Eq. (38). The sum of the gravitational force and the pressure force (or buoyant force, which acts parallel to the external force but in the opposite direction) is written as: ! ρg ~ ~ ~g , Fgra + F pre = m p 1 − (A.1) ρp where ~g is the gravitational acceleration vector. The drag force is given by Sommerfeld (2003)   π F~dra = C D ρg d2p u~g − u~p u~g − u~p 8 (A.2) with  i 24 h  0.687   ; Re p ≤ 1000 1 + 0.15 Re  p  Re p CD =      0.44 ; Re p > 1000 (A.3) The lift force is based on the analytical result of Saffman (1965). For a high particle Reynolds number, this force is given by Crowe et al. (1998): s     ρg µg 2 ~ ~ g × u~g − u~p f Re p , ReS , F sa f = 1.615 d p (A.4) ω ω~g ~ g denotes the angular velocity vector of the fluid phase, defined as: where ω 1 ~ g = ∇ × ~ug , ω 2 (A.5) and the empirical factor is given by Sommerfeld (2003) and Götz (2006) s s   !     Re Re ReS p  S    exp − 1 − 0.3314 + 0.3314      2Re 10 2Re  p p f Re p , ReS =   r    ReS     0.0524 2 ; Re p ≤ 40 (A.6) ; Re p > 40 where ReS represents the Reynolds number of the shear flow., defined as: ReS = ~g ρg d2p ω µg . (A.7) An additional force (called the Magnus force) is exerted on the particles when the particle rotation rate differs from that of the surrounding fluid, given by Crowe et al. (1998)   ~ rel × u~g − u~p Ω π F~mag = Cmag ρg d2p u~g − u~p (A.8) 8 ~ rel Ω ~ rel represents the relative angular velocity between the fluid and solid phases Ω ~ rel = ω ~p − ω ~ g , and the lift where Ω coefficient Cmag is expressed according to Götz (2006) and Alobaid (2013)     0.4G ; G ≤ 1 (A.9) Cmag =    0.4 ; G>1 32 Powder Technology 264 (2014) 343–364 with ~ rel dp Ω G= 2 u~g − u~p . (A.10) The resulted torque acting on a particle given by Eq. (37) is calculated as a sum of two contributions, namely the ~ rot and the moment of force due to shortmoment of force due to the viscous interaction of a particle with the fluid M t : ~ con range tangential contact force M t ~p = M ~ rot + M ~ con M . (A.11) ~ rot is given by Crowe et al. (1998) and Sommerfeld The torque acting on a rotating particle due to a local fluid rotation M (2003): ~ rot M ρg dp = CR 2 2 !5 ~ rel ~ rel Ω Ω (A.12) with a rotational coefficient CR based on the numerical simulations of Dennis et al. (1980) and experimental data of Sawatzki (1970) for high particle Reynolds numbers found to be:  64π        ReR CR =   12.9 128.4      Re0.5 + ReR R ; ReR ≤ 32 (A.13) ; 32 < ReR < 1000 where ReR is the rotating Reynolds number ReR defined as: ReR = ~ rel ρg d2p Ω µg . (A.14) The second term on the right-hand side of Eq. (A.11) is given by: t ~ con M = k=N X k=1 t ri F~con,k × ~ni j (A.15) where N represents the total number of collision partners, which are in contact with this particle. 33 Powder Technology 264 (2014) 343–364 Appendix B. Euler-Euler Model (Two-Fluid Model) Gas Phase Mass flow rates Density Dynamic viscosity Particle Phase Density Diameter Granular viscosity Granular bulk viscosity Granular temperature Solids pressure Radial distribution Maximum packing limit Drag models Restitution coefficients Specularity coefficients Solution Methods & Solver Settings Pressure–velocity coupling Spatial discretization Temporal discretization Time step 0.005 & 0.006 kg/s 1.225 kg/m3 1.79 × 10−5 kg/(m · s) 2500 kg/m3 0.0025 m (monodisperse) Syamlal et al. (1993) Lun et al. (1984) Algebraic & Partial differential approach Lun et al. (1984) Lun et al. (1984) 0.63 Syamlal and O’Brien (1989) Gidaspow et al. (1992) Wen and Yu (1966) 0.6, 0.7, 0.8 & 0.9 0.0, 0.1, 0.2, 0.5, 0.7 & 1.0 Phase Coupled SIMPLE First-order upwind First-order implicit 2 × 10−4 s Table B.2: Simulation conditions and input parameters for the two-fluid model. 34 Powder Technology 264 (2014) 343–364 Appendix C. Euler-Lagrange Model (Discrete Element Method) Gas Phase Mass flow rates Density Dynamic viscosity Particle Phase Density Diameter Total number of particles (P – P) RC in tangential direction (P – P) RC in normal direction (P – W) RC in tangential direction (P – W) RC in normal direction (P – P) and (P – W) friction coefficients Damping parameter in normal direction Damping parameter in tangential direction (P – P) and (P – W) stiffness coefficients Particle sub–time step Solution Methods & Solver Settings Under-relaxation factor → velocity Under-relaxation factor → pressure Total time Solver 0.005 & 0.006 kg/s 1.225 kg/m3 1.79 × 10−5 kg/(m · s) 2500 kg/m3 0.0025 m (monodisperse) 36,500 0.33 Kharaz et al. (1999, 2001) 0.97 Kharaz et al. (1999, 2001) 0.33 Kharaz et al. (1999, 2001) 0.97 Kharaz et al. (1999, 2001) 0.10 Kharaz et al. (1999, 2001) calculated using Eq. (45) calculated using Eq. (50) 4.1 × 105 N/m Di Renzo and Di Maio (2004) 2 × 10−4 s 0.3 0.2 500 ms Multi–Grid/red–black–SOR Table C.3: Simulation conditions and input parameters for the discrete element method; P=Particle , W=Wall and RC=restitution coefficient. 35 Powder Technology 264 (2014) 343–364 References Alobaid, F., 2013. 3D Modeling and simulation of reactive fluidized beds for conversion of biomass with discrete element method. Ph.D. thesis. Technische Universität Darmstadt, Germany. Alobaid, F., Baraki, N., Epple, B., 2014. Investigation into improving the efficiency and accuracy of CFD/DEM simulations. Particuology . Alobaid, F., Epple, B., 2013. Improvement, validation and application of CFD/DEM model to dense gas–solid flow in a fluidized bed. Particuology 11, 514–526. Alobaid, F., Ströhle, J., Epple, B., 2013. Extended CFD/DEM model for the simulation of circulating fluidized bed. Advanced Powder Technology 24, 403–415. Armenio, V., Fiorotto, V., 2001. The importance of the forces acting on particles in turbulent flows. Physics of Fluids 13, 2437–2440. Bokkers, G.A., Van Sint Annaland, M., Kuipers, J.A.M., 2004. Mixing and segregation in a bidisperse gas–solid fluidized bed: a numerical and experimental study. Powder Technology 140, 176–186. Chiesa, M., Mathiesen, V., Melheim, J.A., Halvorsen, B., 2005. Numerical simulation of particulate flow by the Eulerian–Lagrangian and the Eulerian–Eulerian approach with application to a fluidized bed. Computers and Chemical Engineering 29, 291–304. Crowe, C.T., Sommerfeld, M., Tsuji, Y., 1998. Multiphase flows with droplets and particles. CRC Press. Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies. Geotechnique 29, 47–65. Dalla Valle, J.M., 1948. Micromeritics. Pitman, London. Deen, N.G., Van Sint Annaland, M., Van der Hoef, M.A., Kuipers, J.A.M., 2007. Review of discrete particle modeling of fluidized beds. Chemical Engineering Science 62, 28–44. Dennis, S.C.R., Singh, S.N., Ingham, D.B., 1980. The steady flow due to a rotating sphere at low and moderate Reynolds numbers. Journal of Fluid Mechanics 101, 257–279. Di Renzo, A., Di Maio, F.P., 2004. Comparison of contact-force models for the simulation of collisions in DEM-based granular flow codes. Chemical Engineering Science 59, 525–541. Du, W., Bao, X., Xu, J., Wei, W., 2006a. Computational fluid dynamics (CFD) modeling of spouted bed: Assessment of drag coefficient correlations. Chemical Engineering Science 61, 1401–1420. Du, W., Bao, X., Xu, J., Wei, W., 2006b. Computational fluid dynamics (CFD) modeling of spouted bed: Influence of frictional stress, maximum packing limit and coefficient of restitution of particles. Chemical Engineering Science 61, 4558–4570. Duarte, C.R., Olazar, M., Murata, V.V., Barrozo, M.A.S., 2009. Numerical simulation and experimental study of fluid–particle flows in a spouted bed. Powder Technology 188, 195–205. Elghobashi, S.E., Abou-Arab, T..W., 1983. A Two-Equation Turbulence Model for Two-Phase Flows. Phys. Fluids 26, 931–938. Ergun, S., 1952. Fluid flow through packed columns. Chemical Engineering Progress 48, 89–94. Geldart, D., 1973. Types of Gas Fluidization. Powder Technology 7, 285–292. Gera, D., Gautam, M., Tsuji, Y., Kawaguchi, T., Tanaka, T., 1998. Computer simulation of bubbles in large-particle fluidized beds. Powder Technology 98, 38–47. Gidaspow, D., Bezburuah, R., Ding, J., 1992. Hydrodynamics of circulating fluidized beds, kinetic theory approach in fluidization. Proceedings of the 7th Engineering Foundation Conference on Fluidization , 75–82. Goldschmidt, M.J.V., Beetstra, R., Kuipers, J.A.M., 2004. Hydrodynamic modelling of dense gas-fluidized beds: comparison and validation of 3D discrete particle and continuum models. Powder Technology 142, 23–47. Götz, S., 2006. Gekoppelte CFD-DEM-Simulation blasenbildender Wirbelschichten. Shaker Verlag. Gryczka, O., Heinrich, S., Deen, N.G., Annaland, M.v.S., Kuipers, J.A.M., Jacob, M., Mörl, L., 2009. Characterization and CFD–modeling of the hydrodynamics of a prismatic spouted bed apparatus. Chemical Engineering Science 64, 3352–3375. Hertz, H., 1882. Über die Berührung fester elastischer Körper. Journal für die Reine und Angewandte Mathematik 92, 156–171. Hill, R.J., Koch, D.L., Ladd, J.C., 2001. Moderate-Reynolds-numbers flows in ordered and random arrays of spheres. Journal of Fluid Mechanics 448, 243–278. Hjelmfelt, A.T., Mockros, L.F., 1966. Motion of discrete particles in a turbulent fluid. Applied Scientific Research 16, 149–161. Van der Hoef, M.A., van Sint Annaland, M., Deen, N.G., Kuipers, J.A.M., 2008. Numerical simulation of dense gas-solid fluidized beds: a multiscale modeling strategy. Annual Reviews Fluid Mechanics 40, 47–70. Huilin, L., Yurong, H., Wentie, L., Ding, J., Gidaspow, D., Bouillard, J., 2004. Computer simulations of gas–solid flow in spouted beds using kinetic-frictional stress model of granular flow. Chemical Engineering Science 59, 865–878. Johnson, P.C., Jackson, R., 1987. Frictional-collisional constitutive relations for granular materials, with application to plane shearing. Journal of Fluid Mechanics 176, 67–93. Kharaz, A., Gorham, D., Salman, A., 1999. Accurate measurement of particle impact parameters. Measurement Science and Technology 10, 619–647. Kharaz, A., Gorham, D., Salman, A., 2001. An experimental study of the elastic rebound of spheres. Powder Technology 120, 281–291. Lan, X., Xu, C., Gao, J., Al-Dahhan, M., 2012. Influence of solid–phase wall boundary condition on CFD simulation of spouted beds. Chemical Engineering Science 69, 419–430. Link, J., 2006. Development and validation of a discrete particle model of a spout-fluid bed granulator. Ph.D. thesis. University of Twente, The Netherlands. Lun, C.K., Savage, S.B., Jeffrey, D.J., Chepurniy, N., 1984. Kinetic theories for granular flow: inelastic particles in Couette flow and slightly inelastic particles in general flow field. Journal of Fluid Mechanics 140, 223–256. Pei, P., Zhang, K., Wen, D., 2011. Comparative analysis of CFD models for jetting fluidized beds: The effect of inter-phase drag force. Powder Technology 221, 114–122. 36 Powder Technology 264 (2014) 343–364 Saffman, P.G., 1965. The lift on a small sphere in a slow shear flow. Journal of Fluid Mechanics 22, 385–400. Sawatzki, O., 1970. Das Strömungsfeld um eine rotierende Kugel. Acta Mechanica 9, 159–214. Schaeffer, D., 1987. Instability in the Evolution Equations Describing Incompressible Granular Flow. Journal of Differential Equations 66, 19–50. Sommerfeld, M., 2003. Analysis of collision effects for turbulent gas–particle flow in a horizontal channel: Part I. Particle transport. International Journal of Multiphase Flow 29, 675–699. Syamlal, M., O’Brien, T.J., 1989. Computer simulation of bubbles in a fluidized bed. American Institute of Chemical Engineers Symposium Series 85, 22–31. Syamlal, M., Rogers, W., O’Brien, T.J., 1993. MFIX Documentation: Theory Guide. National Technical Information Service, Springfield, VA 1. ANSYSr , 2012. Ansys fluent 12.0 theory guide . Tsuji, T. Yabumoto, K., Tanaka, T., 2008. Spontaneous structures in three-dimensional bubbling gas-fluidized bed byparallel DEM-CFD coupling simulation. Powder Technology 184, 132–140. Tsuji, Y., Kawaguchi, T., Tanaka, T., 1992. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. Powder Technology 71, 239–250. Tsuji, Y., Kawaguchi, T., Tanaka, T., 1993. Discrete particle simulation of 2–dimensional fluidized–bed. Powder Technology 1, 79–87. Wen, C.Y., Yu, Y.H., 1966. Mechanics of fluidization. Chemical Engineering Progress Symposium Series , 100–111. Zhu, H.P., Zhou, Z.Y., Yang, R.Y., Yu, A.B., 2008. Discrete particle simulation of particulate systems: a review of major applications and findings. Chemical Engineering Science 63, 5728–5770. 37