A Comparative Study of Different CFD-codes for Fluidized Beds

1. Introduction

Fluidized beds are one of the most widely used chemical reactors in many industry areas. It has been used mostly for gas-solid reactions as it provides some advantages including its easiness of handling and transport of solids due to fluid-like behavior, intensive mixing of solids leads to a consistent temperature distribution and favorable conditions for heat and mass transfer. [1]. Despite its advantages, it has also some drawbacks which heavily relies upon then complexity of hydrodynamics behavior and the difficulty in scaling up the reactor [1,2]. Historically, empirical methodologies derived from laboratory-scale investigations of fluidized-bed reactors have been instrumental in calculating key operational parameters, such as minimum fluidization velocity, bubble dynamics, particle attrition rates, and the axial distribution of solid mass fractions [3].

Numerical simulations offer a shortcut to circumvent the extensive process involved in designing and executing experiments, enabling rapid assessment of relevant fields at the industrial scale [4]. One of the numerical methods widely used is so-called Computational fluid dynamics (CFD). Several CFD models exist for simulating gas-solid flows, including the two-fluid model, discrete-particle model, hybrid model, and direct numerical simulations. The two-fluid method, particularly its most utilized form, the Eulerian-Eulerian (E-E) approach, considers both phases as interpenetrating continua. Closure relations based on the kinetic theory of granular flows (KTGF) help describe the solid phase with gas-like properties, considering the particles’ motion influenced by interactions within the gas phase and with other particles or walls [5]. In the discrete-particle approach, the fluid phase is treated as a continuum, while the solid phase is considered discrete, tracking individual particles within the computational domain. This method, notably utilized in the Eulerian-Lagrangian (E-L) framework, is recognized for its unique particle collision detection, which can be either stochastic or deterministic, facilitating detailed analysis of particle dynamics [6]. In the hybrid methods, the two-fluid and discrete-particle are integrated to leverage both models’ benefits. This sub-domain method simplifies the computational domain based on an equilibrium index, optimizing simulation efficiency and accuracy [7]. In Direct Numerical Simulation (DNS), it solves directly the Navier-Stokes equation aiming to capture the full range of scales and dynamics by using high-resolution grids and accurate numerical schemes. However, DNS is also limited by its high computational cost and low Reynolds numbers [8].

Over the past 20 years, numerous studies on gas-solid CFD models have been conducted using commercial CFD software, e.g., ANSYS FLUENT [9], ANSYS CFX [10], OpenFOAM [11], MFiX [12], BARRACUDA Virtual Reactor [13], CFDEM coupling [14], COMSOL Multiphysics [15], STAR-CCM [16] and EDEM [17]. Given the fluidized bed’s complex hydrodynamics, researchers focus on examining various parameters that influence gas-solid flow dynamics. Essential factors to consider in numerical models include superficial gas velocity, drag, specularity and restitution coefficients, and frictional viscosity, as these significantly affect the behavior of fluidized beds [18].

Taghipour et al. [19] compared hydrodynamics in a fluidized bed through experiments and FLUENT 6.0 simulations. The study used a pseudo-2D experimental setup and a 2D Eulerian-Eulerian model incorporating Gidaspow [5], Syamlal-O`Brien [20], and Wen-Yu drag models. Results showed the bed reaching steady state in 3 seconds, with consistent bed expansion and bubble formation. While the drag models closely matched experimental observations, discrepancies in voidage and bed expansion highlighted the need for further CFD model validation.

Cornelissen et al. [21] investigated the impact of mesh size, time step, and convergence criteria on the accuracy of a liquid-solid fluidized bed model in FLUENT 6.1.22, benchmarking their findings against experimental data. They determined that a convergence criterion of ${10}^{-4}$ led to faster, convergent results compared to the standard ${10}^{-3}$. Additionally, they noted that Courant numbers within the range of $0.03$ to $0.3$ produced reliable outcomes regardless of mesh size and time step.

Li et al. [22] explored how solid phase boundary conditions affect simulations by adjusting specularity and particle-wall restitution coefficients in both 2D and 3D Eulerian-Eulerian models using MFIX software, comparing results to a pseudo 2D experiment. They discovered that the restitution coefficient’s impact was minimal compared to the specularity coefficient, with specularity coefficients under 0.05 leading to notable deviations in solid flow near walls and altering circulation patterns. While both 2D and 3D simulations captured the general flow behavior, 3D simulations aligned more closely with experimental data, emphasizing the significance of wall effects in pseudo 2D experiments for CFD model validation.

Herzog et al. [23] conducted a comparative analysis using MFIX, OpenFOAM versions 1.6/2.0, and FLUENT 6.3 against experimental findings by Taghipour, focusing on a 2D Eulerian-Eulerian using Gidaspow and Syamlal O`Brien drag models. The study revealed that FLUENT and MFIX closely matched experimental data for both drag models, whereas OpenFOAM’s predictions were not satisfactory, suggesting OpenFOAM requires further refinement and validation for fluidization prediction.

Shi et al. [18] carried out the research to align 2D and 3D models with Taghipour et al. [19] experimental findings, utilizing Eulerian-Eulerian frameworks in FLUENT 16.2. Their investigation spanned dimensions, flow regimes, boundary conditions, mesh granularity, and drag models. Optimal grid size emerged as 18 times the particle diameter for closest experimental alignment, with finer meshes diminishing accuracy. Minor variances were observed within specific specularity and restitution coefficient ranges, deeming them suitable. Turbulence’s impact was minimal within the bed but influenced gas flow above, recommending laminar models for non-reactive, isothermal beds and RANS $k-ϵ$ models otherwise. Syamlal-O`Brien drag models, particularly with Johnson-Jackson boundary conditions and turbulent RANS models, showed superior experimental congruence. They noted 3D models matched experimental data more closely than 2D, cautioning that certain 2D model parameter combinations might mislead due to unrealistic assumptions. Consequently, they advocated for 3D validations, reserving 2D analyses for pre-validated model sensitivity studies.

In the current study, the fluidized bed columns from Taghipour et al. [19] were analyzed in both 2D and 3D simulations using three different CFD software codes: ANSYS Fluent 2023 R3, MFiX 23.1.1, and OpenFOAM 10. Simulations in MFiX and OpenFOAM were realized as part of the MSc thesis of Daun [24]. The current work extends the analysis presented by [24], and compares results further with new simulations realized in ANSYS Fluent. This investigation included the examination of two drag models, Syamlal-O’Brien and Gidaspow, across varied superficial gas velocities. Additionally, the study explored the effects of varying restitution and specularity coefficients. Simulations employed a Eulerian-Eulerian framework underpinned by the kinetic theory of granular flow (KTGF) aiming to highlight and address the subtleties inherent in inter-software comparisons.

2. Problem Description

The observation of several parameters as outlined previously was done by replicating the experimental case from Taghipour et al. [19] as it enables modelling validation by comparing to the experimental results therein. The model’s geometry consists of a rectangular column, measuring 0.28 m in width, 1 m in height, and 0.025 m in depth. The lower section of the column, extending up to 0.4 meters, is densely packed with granular particles with mean particle size of 275 $\mu$m at a packing density of 0.6, illustrated in Figure 1. According to Taghipour et al. [19] the minimum fluidization velocity ($U_{mf}$) is 0.065 m/s, whereas the system exhibited an overall pressure drop of 4400 Pa, a bed expansion ratio of 1.1, and a voidage (gas-phase volume fraction) of 0.5.

3. Computational Model

In this work, the simulations were carried out in the framework of Eulerian-Eulerian (E-E) by utilizing the kinetic theory of granular flow (KTGF). Three different CFD codes, namely, ANSYS Fluent version 2023 R3, MFiX 23.1.1 and OpenFOAM 10 were evaluated along with the implementation of Syamlal-O`Brien and Gidaspow drag model by varying superficial gas velocities up to 0.03, 0.1, 0.38, 0.46, and 0.51 m/s, equivalent to 0.5, 1.5, 6, 7, and 8 times the minimum fluidization velocity ($U_{mf}$). The complete selected simulation model parameters are summarized in Table 1. The governing equation and simulation setting of ANSYS Fluent, MFiX and OpenFOAM for developing the fluidized bed simulation model are explained below.

3.1. Governing Equation

The governing equations used in this simulation of fluidized bed for air and glass beads case can be summarized as follows:

Mass balance equations

$$\frac{\partial }{\partial t}\left({\alpha }_g{\rho }_g\right)+\nabla ·\left({\alpha }_g{\rho }_g{\overrightarrow{u}}_g\right)=0$$

(1)

$$\frac{\partial }{\partial t}\left({\alpha }_s{\rho }_s\right)+\nabla ·\left({\alpha }_s{\rho }_s{\overrightarrow{u}}_s\right)=0$$

(2)

Where ${\alpha }_g$, ${\alpha }_s$, ${\overrightarrow{u}}_g$, ${\overrightarrow{u}}_s$, ${\rho }_g$ and ${\rho }_s$ represent the volume fraction, velocity and density of gas phase and solid phase respectively. The total volume fraction of gas and solid phase is equal to unity:

$${\alpha }_g+{\alpha }_s=1$$

(3)

Momentum balance equations

$$\frac{\partial }{\partial t}\left({\alpha }_g{\rho }_g{\overrightarrow{u}}_g\right)+\nabla ·\left({\alpha }_g{\rho }_g{\overrightarrow{u}}_g{\overrightarrow{u}}_g\right)=-{\alpha }_g\nabla p+\nabla ·{\overline{\overline{\tau }}}_g+{\alpha }_g{\rho }_g\overrightarrow{g}+K_{gs}\left({\overrightarrow{u}}_g-{\overrightarrow{u}}_s\right)$$

(4)

$$\frac{\partial }{\partial t}\left({\alpha }_s{\rho }_s{\overrightarrow{u}}_s\right)+\nabla ·\left({\alpha }_s{\rho }_s{\overrightarrow{u}}_s{\overrightarrow{u}}_s\right)=-{\alpha }_s\nabla p-\nabla p_s+\nabla ·{\overline{\overline{\tau }}}_s+{\alpha }_s{\rho }_s\overrightarrow{g}+K_{gs}\left({\overrightarrow{u}}_g-{\overrightarrow{u}}_s\right)$$

(5)

Where $K_{gs}$ represents the momentum exchange coefficient or interphase drag force coefficient which will be given by following Syamlal-O`Brien and Gidaspow drag models. The terms ${\overline{\overline{\tau }}}_g$ and ${\overline{\overline{\tau }}}_s$ represent the deviatoric stress tensors for the gas and solid phases, respectively. Pressure gradients are accounted for by $\nabla p$ in the gas phase and $\nabla p_s$ in the solid phase. The gravitational force is represented by $\overrightarrow{g}$. The velocity of gas phase and solid phase are given by ${\overrightarrow{u}}_g$ and ${\overrightarrow{u}}_s$.

Syamlal-O`Brien drag model [20]

$$K_{gs}=\frac{3}{4}\frac{{\alpha }_s{\alpha }_g{\rho }_g}{v_{r,s}^2{d}_s}{C}_D\left(\frac{R{e}_s}{u_{r,s}}\right)\left|{\overrightarrow{u}}_s-{\overrightarrow{u}}_g\right|$$

(6)

where

$$C_D={\left(0.63+\frac{4.8}{\sqrt{R{e}_s/u_{r,s}}}\right)}^2$$

(7)

and

$$u_{r,s}=0.5\left(A-0.06R{e}_s+\sqrt{{\left(0.06R{e}_s\right)}^2+0.12R{e}_s(2B-A)+A^2}\right)$$

(8)

where

$$R{e}_s=\frac{{\rho }_g{d}_s\left|{\overrightarrow{u}}_s-{\overrightarrow{u}}_g\right|}{{\mu }_g}$$

(9)

with

A=${\alpha }_g^{4.14}$, B=$0.8{\alpha }_g^{1.28}$ for ${\alpha }_g⩽0.85$, or

A=${\alpha }_g^{4.14}$, B=${\alpha }_g^{2.65}$ for ${\alpha }_g>0.85$

Gidaspow drag model [5]

For ${\alpha }_g>0.8$

$$K_{gs}=\frac{3}{4}{C}_D\frac{{\alpha }_s{\alpha }_g{\rho }_g\left|{\overrightarrow{u}}_s-{\overrightarrow{u}}_g\right|}{d_s}{\alpha }_g^{-2.65}$$

(10)

for ${\alpha }_g⩽0.8$

$$K_{gs}=150\frac{{\alpha }_s^2{\mu }_g}{{\alpha }_g{d}_s^2}+1.75\frac{{\alpha }_s{\rho }_g\left|{\overrightarrow{u}}_s-{\overrightarrow{u}}_g\right|}{d_s}$$

(11)

where

$$C_D=\frac{24}{{\alpha }_{g}R{e}_s}\left(1+0.15{\left({\alpha }_{g}R{e}_s\right)}^{0.687}\right)$$

(12)

Solid phase fluctuating energy equations [28]

$$\frac{3}{2}\left[\frac{\partial }{\partial t}\left({\rho }_s{\alpha }_s{\Theta }_s\right)+\nabla ·\left({\rho }_s{\alpha }_s{\overrightarrow{u}}_s{\Theta }_s\right)\right]=\left(-p_s\overline{I}+{\overline{\tau }}_s\right):\nabla {\overrightarrow{u}}_s+\nabla ·\left(k_{{\Theta }_s}\nabla {\Theta }_s\right)-{\gamma }_{{\Theta }_s}+{\phi }_{gs}$$

(13)

In order to attain closure of the solid stress terms in equation 4, 5 and 13 KTGF was used to derive the constitutive equations provided in in Table A1. KTGF provides values for solid-phase shear and bulk viscosities, as well as solid pressure. This theory evolves from Ogawa et al. [29]. It assumes particles are smooth, with inelastic and instantaneous binary collisions, as detailed by Gidaspow [5]. Drag models from Syamlal and O’Brien [20] and Gidaspow [5] were used in this study to calculate the momentum exchange coefficients as shown in equation 6 to 12. Conversely, the Gidaspow drag model merges the Wen-Yu drag model with the Ergun equation that is pertinent for dense beds, linking drag to pressure drop across porous media.

3.2. Simulation Setup

The OpenFOAM simulations adopted the fluidized bed tutorial case on multiphaseEulerFoam for fluidized Bed, whereas MFiX simulations used the 2D TFM fluid bed preset, and ANSYS Fluent based its setup on a 2D fluidized bed case from the training materials [30]. The computational domains for 2D and 3D models were discretized into 12,200 and 56,000 rectangular cells, respectively, with a grid width of $\Delta x=0.005$ m. Boundary conditions for all codes were set to default. Inlet velocities were varied, 0.03, 0.1, 0.38, 0.46, and 0.51 m/s. For the gas phase at the wall, a no-slip condition for the gas phase was applied at the walls. All three CFD codes facilitate the application of the Johnson and Jackson partial slip boundary condition for particles, enabling nuanced simulation of particle-wall interactions [31]. This formulation, articulated by Johnson and Jackson [31], prescribes boundary conditions for velocity and granular temperature in the solid phase, satisfying the specified requirements as outlined in equation 14 and 15. The default specularity coefficient of 0.2 and restitution coefficient of 0.9 for all codes.

$${\mu }_s\frac{\partial {\overrightarrow{u}}_s}{\partial x}=\frac{\pi K{\rho }_s{\alpha }_s{g}_{0,ss}\sqrt{{\Theta }_s}}{2\sqrt{3}{\alpha }_s^{max}}{\overrightarrow{u}}_s$$

(14)

$${\kappa }_s\frac{\partial K}{\partial x}=-\frac{\pi K{\overrightarrow{u}}_s^2{\rho }_s{\alpha }_s{g}_{0,ss}\sqrt{{\Theta }_s}}{2\sqrt{3}{\alpha }_s^{max}}-\frac{\pi \sqrt{3}{\rho }_s{\alpha }_s{g}_{0,ss}\left(1-e_{ss}^2\right)\sqrt{{\Theta }_s}}{4{\alpha }_s^{max}}{\Theta }_s$$

(15)

Where K and $e_{ss}$ represent the specularity and restitution coefficients respectively, their values are constrained between 0 and 1. A specularity coefficient (K) of 0 denotes perfectly specular collisions, whereas a value of 1 indicates perfectly diffuse collisions. Similarly, a restitution coefficient ($e_{ss}$) of 0 corresponds to perfectly inelastic collisions, and a value of 1 represents perfectly elastic collisions. Motivated by the clear graphical analysis provided by [18], which elucidates the influence of these coefficients on fluidized bed flow behavior, this study set forth to investigate the effects of several restitution (0.9, 0.95, and 0.99) and specularity (0.05, 0.1, and 0.2) coefficient values to evaluate their impact on simulation results. In total, 50 simulation cases were executed, employing specific parameters i.e. CFD codes, superficial gas velocities, drag models, restitution coefficient, specularity coefficient and effect of 2D and 3D geometry as delineated in Table 3, to thoroughly investigate their effects.

Based on the sensitivity analysis conducted by Taghipour et al. [19] regarding time step and spatial discretization, minimal differences were observed between time steps of 0.0005s and 0.001s, and between second and first-order discretizations. Consequently, this study adopted a time step of 0.001s. In the discretization settings, ANSYS Fluent employed a first-order upwind and PRESTO for spatial discretization, alongside a first-order implicit method for temporal discretization. MFiX utilized a high order upwind scheme with a downwind factor for spatial discretization, incorporating the Superbee flux limiter. This approach was complemented by an implicit Euler method for temporal aspects. OpenFOAM’s strategy involved Gaussian discretization, applying flux limiters like "vanLeer" or "limitedLinear" based on the variable, highlighting a targeted approach to discretization.

The numerical settings included the use of the AMG solver coupled with the SIMPLE algorithm in ANSYS Fluent to address the software`s fluid flow challenges. MFiX and OpenFOAM utilized the PBiCGStab solver. MFiX employed a specific algorithm designed for its simulations [32], whereas OpenFOAM used the PIMPLE algorithm, indicating distinct computational strategies adopted by each software to navigate the complexities of multiphase flow simulations. Regarding convergence criteria, ANSYS Fluent was set to default settings for all variables at a tolerance of 0.001. OpenFOAM adopted strict tolerances, with 1e-9 for volume fraction (alpha.*), 1e-8 for pressure head (p_rgh) and velocity (U.*), and for temperature-related fields (Theta.*), with turbulence kinetic energy and its dissipation rate (k|epsilon.*) at 1e-5. MFiX set its default continuity tolerance at 0.001, with tolerances for energy, species, granular energy, and scalar variables at 0.0001. Under-relaxation factors were kept at default settings for ANSYS Fluent, OpenFOAM, and MFiX, reflecting a standard approach to ensuring numerical stability within each software framework.

The simulation duration was set to 12 seconds, with time-averaging from 3 to 12 seconds, as suggested by Herzog et al. [23] and Taghipour et al. [19], to capture a stable and representative state of the system being analyzed.

4. Results

The results, derived from various simulations, are delineated in terms of bed expansion ratio ($H/H_0$), pressure drop, solid volume fraction contour plots, solid velocity at $y=0.2$ m, and voidage at $y=0.2$ m. These results were obtained using three distinct CFD codes and evaluated across two drag models. The bed expansion ratio and pressure drop data encompass all examined superficial velocities: 0.03, 0.1, 0.38, 0.46, and 0.51 m/s, corresponding to 0.5, 1.5, 6, 7, and 8 times $U_{mf}$, respectively. In contrast, the solid volume fraction contour plots specifically showcase three selected superficial velocities: 0.2, 0.38, and 0.46 m/s, at times of 0.4 s, 1 s, and 3 s. Additionally, the analyses of solid velocity and voidage at $y=0.2$ m encapsulate the effects of variations across all two drag models, and the three CFD codes at 0.38 m/s of superficial gas velocity. The simulation outcomes, including bed expansion ratio, pressure drop, and voidage, were evaluated against the experimental results from Taghipour et al. [19], thereby assessing the simulations’ accuracy by using the formula for standard deviation as elucidated in Shi et al. [18].

4.1. Bed Expansion Ratio and Pressure Drop

The expansion of the bed was assessed by charting the time-averaged solid fraction at the bed’s center vertically. A criterion was established where a solid fraction value below 1% indicated the extent of bed expansion. This procedure was applied across various superficial gas velocities to observe how bed expansion varied with changes in inlet velocity. The summary of the results of bed expansion and pressure drop are summarized in Table 1. Figure 2 presents the bed expansions determined using the Syamlal-O`Brien drag model within Fluent, MFiX and OpenFOAM. From this figure, it is evident that all CFD codes forecast minimal bed expansion at lower velocities, with a linear increase in expansion noted for velocities exceeding 0.38 m/s. In 3D scenarios, bed expansion appears reduced compared to analogous 2D situations. Moreover, the figure illustrates that Fluent tends to forecast greater bed expansion compared to MFiX and OpenFOAM. Figure 3 illustrates the bed expansions observed when employing the Gidaspow drag model in MFiX and OpenFOAM. From this figure, it is apparent that three software packages register bed expansion at all velocities above the minimum fluidization velocity. The expansion appears to increase linearly in relation to the superficial gas velocity. The 3D simulations show a notably lesser degree of bed expansion compared to their 2D counterparts. A comparison with Figure 2 reveals that the expansions resulting from the Gidaspow model from OpenFOAM and MFiX exceed those observed with the Syamlal-O`Brien model. Surprisingly, Fluent forecasts lower bed expansions than OpenFOAM and MFiX when utilizing the Gidaspow model. Additionally, a contraction in the bed at the lowest gas velocity is noted in Fluent and MFiX with both drag models, a phenomenon not observed in OpenFOAM simulations.

The pressure drop across the bed was assessed at the outlet of the fluidized bed reactor for each tested superficial gas velocity. Figure 4 and Figure 5 illustrate a pattern where the pressure drop initially spikes to a peak value before gradually diminishing as the superficial gas velocity increases. It is observed that with the Gidaspow model, the peak pressure drop occurs at lower superficial gas velocities compared to the Syamlal-O`Brien model. Additionally, both figures indicate that MFiX and Fluent estimates a greater pressure drop across all tested superficial gas velocities than OpenFOAM, in both 2D and 3D simulations. A comparison between 2D and 3D results shows that Fluent and OpenFOAM reports a smaller pressure drop in 3D compared to 2D, whereas MFiX exhibits an increase in pressure drop from 2D to 3D simulations. It is also worth mentioning, specifically for Fluent, that the pressure drop at superficial velocities lower than the fluidization velocity tends to predict a very high pressure drop. This observation aligns with the findings from Taghipour et al. [19], El Ajouri et al. [33], suggesting that the primary cause could be attributed to the solids not being fully fluidized at these lower velocities, thereby being predominantly influenced by interparticle frictional forces. Such forces are not adequately predicted by the multifluid model used for simulating gas–solid phases in Fluent, which may not incorporate these frictional interactions effectively.

4.2. Solid Velocity and Voidage

Figure 6 displays the time-averaged solid particle velocity in the y-direction at a height of 0.2 meters, with a superficial gas velocity of 0.38 m/s, utilizing both Syamlal-O`Brien drag model in Fluent, MFiX and OpenFOAM. The observed trend across all data lines reveals a downward velocity near the walls transitioning to a parabolic upward velocity near the bed’s center. Notably, Fluent predicts significantly higher velocities than MFiX, particularly near the walls and the reactor’s center. Furthermore, for 3D simulations, velocities were averaged across the z-direction cells to align with the 2D results, offering a comprehensive view of the domain. While MFiX’s 3D outcomes closely resemble its 2D counterparts, OpenFOAM and Fluent markedly diverge, maintaining a relatively constant velocity across the bed as opposed to the parabolic profile observed in 2D simulations. In contrast to the Syamlal-O`Brien drag model, Figure 7 illustrates the time-averaged solid particle velocity using the Gidaspow drag model, processed similarly to Figure 6. The outcomes indicate that Fluent and MFiX yield nearly identical predictions, with velocities marginally higher at the bed’s center but displaying a somewhat flattened profile across both 3D and 2D simulations. Conversely, OpenFOAM’s results diverge slightly, showing higher velocities near the walls and at the center of the bed in both 3D and 2D scenarios. Notably, the 3D results in OpenFOAM exhibit a slightly flatter velocity profile compared to their 2D counterparts, suggesting a distinct behavior in spatial velocity distribution when using the Gidaspow model.

Further results, as shown in Figure 8, compare the time-averaged voidage for 2D and 3D simulations at the same bed height and gas velocity, employing the Syamlal-O’Brien drag model. This figure indicates that Fluent displays the highest voidage in 3D simulations compared to others, presenting an asymmetrical appearance relative to the rest. In OpenFOAM, the 3D simulations show a maximum voidage nearly identical to their 2D counterparts. Conversely, MFiX’s 3D simulations exhibit a zone of low voidage extending from the wall, which sharply increases to a peak significantly greater than that observed in 2D simulations at the center of the bed. On the other hand, Figure 9 showcases the time-averaged voidage for 2D and 3D simulations at the same bed height and gas velocity, utilizing the Gidaspow drag model. From this figure, it is evident that MFiX’s 3D simulations replicate the behavior observed in the Syamlal-O’Brien case, with voidage sharply increasing at the center of the bed in 3D simulations. For the Fluent simulations, both 2D and 3D present almost identical patterns, though the 3D simulations exhibit slightly higher voidage at the left wall. In OpenFOAM, the 3D results are slightly higher than 2D, maintaining a nearly similar pattern.

4.3. Solid Volume Fraction

The evolution of the fluidized beds and time averaged solid volume fractions with three chosen superficial gas velocities equal to 0.2, 0.38, and 0.46 m/s as the results from three different CFD codes with Syamlal-O`Brien drag model as well as Gidaspow drag model are depicted in Figure 10 and Figure 11. Initially, the bed height rose with the formation of bubbles until it stabilized at a constant steady-state height. The solid flow patterns inside the bed are relatively symmetrical for all cases before t < 3 s. The chaotic and transient bubble formation process at 3 seconds indicates that the statistical steady-state condition was achieved. Therefore, it confirmed the previous findings that a simulation duration of 3 seconds is sufficient to attain statistical steady-state behavior across all drag function models [19,23].

At the beginning of simulation, Fluent and MFiX display comparable patterns, with OpenFOAM showing slight deviations. As the flow progresses, the patterns of flow and bubble formation for both Syamlal-O’Brien and Gidaspow models remain similar in Fluent and MFiX but differ in OpenFOAM, highlighting that the numerical results of MFiX and Fluent are quite aligned, as also confirmed by Herzog et al. [23]. A closer examination reveals good agreement in the time-averaged solid volume fraction among all codes. All cases show a lower voidage observed near the wall and at the base of the bed, whereas an increased voidage is noted in the bed’s central area. This pattern of solid volume fraction distribution indicates a higher formation of bubbles in the central region as opposed to the bottom. Specifically, in the Syamlal-O`Brien model, Fluent exhibits a higher voidage and greater bed expansion compared to MFiX and OpenFOAM. Conversely, in the Gidaspow model, OpenFOAM demonstrates larger bed expansion and more voidage. The case studies with three superficial velocities further reveal that MFiX tends to display a more distinct phase border than OpenFOAM and Fluent, where phase interfaces blend more smoothly. This distinction appears to lead to the formation of more, but smaller, bubbles in OpenFOAM, whereas MFiX features fewer but significantly larger bubbles, highlighting the intricate dynamics of fluidized beds as captured by different simulation tools and drag models.

5. Discussion

Fluidized beds continue to be indispensable in the industrial sector, a fact underscored by the extensive volume of research utilizing Computational Fluid Dynamics (CFD) codes for their study, as elaborately discussed by Alobaid et al. [4]. This significant focus on fluidized beds through CFD simulations brings to the forefront the necessity of determining the most suitable CFD software for accurately capturing the complex hydrodynamics within these systems. Among the available CFD tools, ANSYS Fluent and MFiX have been recognized for their maturity and reliability in modeling fluidized bed dynamics. In contrast, OpenFOAM, despite its potential, has faced skepticism regarding its accuracy and applicability for fluidized bed simulations, as highlighted by Herzog et al. (2012). In this context, our study undertook on a comparative evaluation of these three CFD codes by simulating a standard fluidized bed configuration and comparing the outcomes with the experimental findings reported by Taghipour et al. [19]. In addition to that, the examination of each software’s numerical predictions of voidage fraction by employing the standard deviation formula as a measure of accuracy was also presented.

5.1. Comparative study of different codes

Figure 2 and Figure 3 showed a comparison between simulation results and experimental data for Syamlal-O`Brien and Gidaspow drag model, respectively. It is crucial to note that a challenge with the experimental data lies in the absence of a detailed description of how bed expansion measurements were conducted, limiting the direct comparability to the current simulation findings. Despite this, a general consistency in trends between simulated and experimental results is observable. For superficial velocities below 0.38 m/s using the Syamlal-O’Brien drag model, in the case of using Syamlal-O`Brien drag model, ANSYS Fluent tends to overestimate bed expansion, whereas MFiX and Fluent show a slight underestimation. However, starting from a superficial velocity of 0.38 m/s, all three software overestimate bed expansion. It’s noteworthy that results from 3D simulations presented lower bed expansion compared to their 2D counterparts, underscoring the significance of incorporating friction along the front and back walls, regardless of the drag model applied. While MFiX and OpenFOAM align closely with experimental outcomes, ANSYS Fluent consistently overestimates bed expansion, highlighting discrepancies across different CFD tools in predicting fluidized bed dynamics. This can also be clearly observed from the solid phase distribution pattern depicted in Figure 10, which shows that ANSYS Fluent, using the Syamlal-O’Brien drag model, exhibits higher bed expansion and significantly lower particle occupancy, or higher voidage, appearing much more fluidized compared to other software utilizing the same drag model. Conversely, with the Gidaspow drag model, the tendency of simulation results to overestimate bed expansion at and above 0.38 m/s persists. However, in this scenario, ANSYS Fluent typically underpredicts bed expansion, followed closely by MFiX, while OpenFOAM tends to slightly overpredict. Nonetheless, starting from a superficial velocity of 0.38 m/s, all three software solutions overestimate the bed expansion. Despite these tendencies, it is observed that findings from the 3D simulations align well with the experimental data.

The comparison of simulation outcomes with experimental data on pressure drop, employing the Syamlal-O’Brien and Gidaspow models, is illustrated in Figure 4 and Figure 5, respectively. Experimentally, the pressure drop is observed to sharply increase until reaching the minimum fluidization velocity, thereafter ascending more gradually with the rise in superficial gas velocity. Contrarily, numerical simulations diverge from this pattern, showing a sharp increase in pressure drop around the minimum fluidization velocity, followed by a gradual decrease as superficial gas velocity increases.

Notably, previous observations from OpenFOAM, as reported by Herzog et al. [23], indicated a significant overestimation of pressure drop. However, current results demonstrate that OpenFOAM now mirrors the trend seen in MFiX, with a pressure drop approximately 500 Pa lower. The introduction of a third dimension reveals differing impacts on pressure drop predictions between MFiX and OpenFOAM across both drag models. MFiX predicts a higher, and OpenFOAM a lower, pressure drop than their 2D simulations, deviating further from the experimental data. Conversely, with ANSYS Fluent, a distinct discrepancy is noted at a superficial velocity of 0.03 m/s across both drag models as previously explained in previous section. Across the board, ANSYS Fluent tends to overestimate pressure drops initially, then adjusts to underestimations. It is evident that no simulation perfectly aligns with the experimental data, with all predicted void fractions being overly large. Despite this, the outcomes could still be deemed in reasonable concordance with experimental observations. However, ANSYS Fluent gives a better fit to the experimental data at superficial velocities of 0.38 m/s and above. Furthermore, the introduction of the third dimension in ANSYS Fluent even results in a much better fit to the experimental result compare to MFiX and OpenFOAM. Notably, while the Syamlal-O`Brien and Gidaspow models yield broadly similar results, the Syamlal-O`Brien model slightly edges closer to the experimental data than Gidaspow in terms of pressure drop. It also align with the resport from Aboudaoud et al. [34] and Liu and Hinrichsen [35] which implied that Syamlal-O`Brien model is more suitable for pressure drop prediction compared to Gidaspow.

5.2. Simulation Error

To assess the accuracy of the numerical predictions of voidage fraction against experimental data from Taghipour et al. [19], the standard deviation profiles for both 2D and 3D simulations were calculated and are presented in Figure 12. The formula for standard deviation is given as follows:

$${\sigma }_{\mathrm{est}}\%=\sqrt{\frac{\sum {\left(Y_{\mathrm{simulation}}-Y_{\mathrm{experiment}}\right)}^2}{N}}\times 100$$

(16)

where ${\sigma }_{\mathrm{est}}$ represents the estimated standard error, $Y_{\mathrm{experiment}}$ and $Y_{\mathrm{simulation}}$ are the experimental and predicted points of the voidage fraction along the x-coordinate at y = 0.2 m, respectively. The standard deviation serves as a reliable measure of error compared to other metrics since this method is derived from Linear Least Squares Regression commonly used in curve fitting. The standard deviation for 3D and 2D models falls within the range of 5.06 – 11.81 %. The findings suggest that neither 3D nor 2D simulations inherently guarantee closer agreement with experimental data. Interestingly, certain 2D simulations exhibit better alignment with the experimental data than their 3D counterparts, indicating that with appropriate parameter selection, 2D simulations can also be a viable option, as supported by Shi et al. [18].

By examining both the qualitative aspects presented in Figure 8 and Figure 9 and the quantitative data outlined in Figure 12, one can discern the impact of varying software and drag models on simulation outcomes. The results demonstrate that MFiX generally yields lower simulation discrepancies, followed by ANSYS Fluent and OpenFOAM. Notably, in the 3D context, MFiX and ANSYS Fluent exhibit similar levels of simulation error. Regarding the choice of drag model, the Syamlal-O`Brien model tends to produce lower simulation errors when used with MFiX and OpenFOAM, whereas the Gidaspow model appears to minimize simulation errors more effectively when applied in ANSYS Fluent. Consequently, while OpenFOAM holds potential for various applications, its performance in fluidized bed simulations appears to lag slightly behind that of MFiX and ANSYS Fluent, highlighting the importance of selecting the appropriate software and drag model to achieve the best simulation accuracy. This observation is consistent with the most significant differences in flow evolution between OpenFOAM and ANSYS Fluent and MFiX, as illustrated in Figure 10 and Figure 11. A previously discussed critical distinction lies in the phase interfaces within the simulations. OpenFOAM’s tendency to generate numerous and diffuse bubbles suggests that gas flow actively propels particles upwards across almost the entire domain. This phenomenon is visible in time-averaged plots from OpenFOAM, where a narrow region of particle accumulation along the wall is the only area where particles can descend. In contrast, MFiX and ANSYS Fluent simulations depict a different behavior, where smaller bubbles at the bed’s bottom coalesce into relatively larger and fewer distinct bubbles that ascend through the bed’s center.

5.3. The Effect of Restitution and Specularity Coefficient

The restitution coefficient is indicative of the collision behavior between particles in the fluidized bed, while the specularity coefficient represents the nature of collisions between particles and the walls. In this study, restitution coefficient values of 0.9, 0.95, and 0.99 were explored. According to the results presented in Table 1, a higher restitution coefficient leads to increased bed expansion in the Syamlal-O’Brien drag model, although it does not significantly affect the pressure drop. This trend is consistent with the Gidaspow drag model, where increasing the restitution coefficient also increases bed expansion, albeit with a slight decrease (0.013) between 0.9 and 0.95, followed by an increase at 0.99.

Analyzing the simulation results based on voidage measurement at y = 0.2 m at a superficial velocity of 0.38 m/s, it was observed that with both drag models, as the restitution coefficient increases (Syamlal-O’Brien: 10.03, 10.90, 13.57%; Gidaspow: 7.84, 8.40, 11.56%), the simulation error also increases this aligns with the result reported by Shi et al. [18]. This suggests that a restitution coefficient of 0.9 most accurately represents the conditions.

Conversely, regarding the impact of the specularity coefficient, which was varied between 0.05, 0.1, and 0.2, the results in the Syamlal-O’Brien drag model indicate that an increase in the specularity coefficient does not significantly alter bed expansion. However, in the Gidaspow drag model, a slight reduction in bed expansion was observed. The simulation error analysis within the Syamlal-O’Brien model correlates with bed expansion outcomes, showing minimal change across specularity coefficients (Syamlal-O’Brien: 10.51, 10.31, 10.03), suggesting that the influence of specularity coefficient variations on voidage predictions can be neglected, as also noted by Shi et al. [18]. Nonetheless, employing the Gidaspow model (Gidaspow: 7.45, 9.24, 7.84) reveals slight differences in simulation error values, with a notable increase at 0.1, indicating a sensitivity to specularity coefficient variations within the Gidaspow drag model.

6. Conclusions

In the current study, the gas-solid flow within 2D and 3D fluidized bed columns, replicated from Taghipour et al. [19], was simulated using the Eulerian-Eulerian framework coupled with the Kinetic Theory of Granular Flow (KTGF). This investigation aimed to assess the capabilities of three prominent CFD software packages—ANSYS Fluent, MFiX, and OpenFOAM—in modeling fluidized bed dynamics. The study specifically focused on the application of the Syamlal-O’Brien and Gidaspow drag models across varying superficial velocities, restitution coefficients, and specularity coefficients to provide a comprehensive evaluation of these tools. The key findings of this research are summarized as follows:

• The assessment of ANSYS Fluent, MFiX, and OpenFOAM showed they could achieve reasonable agreement with experimental data across both drag models, offering insights into fluidized bed hydrodynamics. This comparison underscored specific strengths and weaknesses, revealing MFiX and ANSYS Fluent as more consistent and accurate in simulating fluidized bed dynamics, thereby emerging as preferable choices. OpenFOAM has shown maturity in predicting fluidized bed phenomena, yet it necessitates further validation and exploration to fully grasp the impact of various numerical settings and boundary conditions on its outcomes. This approach underscores the importance of fine-tuning and accurately assessing OpenFOAM’s capabilities for enhanced reliability in fluidized bed simulations.
• The analysis of the Syamlal-O’Brien and Gidaspow drag models revealed significant differences in their predictions, particularly in terms of bed expansion and pressure drop. The Syamlal-O’Brien model, when used with MFiX and OpenFOAM, typically resulted in lower simulation errors compared to the Gidaspow model. However, ANSYS Fluent showed a preference for the Gidaspow model, indicating that the choice of drag model can significantly influence simulation outcomes and should be carefully considered based on the CFD software being used.
• The transition from 2D to 3D simulations introduced notable variations in bed expansion and pressure drop predictions. While 3D simulations offer a more detailed representation of fluidized bed dynamics, they do not always guarantee improved alignment with experimental data. In some cases, 2D simulations provided better fits, suggesting that the dimensionality of the simulation should be chosen based on specific study requirements and the expected flow behavior.
• The study showed that while higher restitution coefficients led to increased bed expansion, indicating their significant impact on particle-particle collisions, variations in the specularity coefficient, which affects particle-wall interactions, had a subtler influence on bed dynamics, particularly noted in the Gidaspow drag model. A restitution coefficient of 0.9 was identified as optimal, balancing accuracy and computational efficiency, whereas the specularity coefficient’s effect, though present, was less pronounced within the examined range. This emphasizes the necessity of carefully selecting these coefficients to accurately model the complex behaviors in fluidized beds.

Appendix A.1

The constitutive equations derived from KGTF is provided in below Table.

Table A1. Constitutive equations used in this study

Gas phase shear stress tensor

$${\overline{\overline{\tau }}}_g={\alpha }_g{\mu }_g\left[\nabla {\overrightarrow{u}}_g+{\left(\nabla {\overrightarrow{u}}_g\right)}^T-\frac{2}{3}\left(\nabla ·{\overrightarrow{u}}_g\right)\overline{\overline{I}}\right]$$ (A1)

Solid phase shear stress tensor

$${\overline{\overline{\tau }}}_s={\alpha }_s{\mu }_s\left[\nabla {\overrightarrow{u}}_s+{\left(\nabla {\overrightarrow{u}}_s\right)}^T+\left(\frac{2}{3}-\frac{{\lambda }_s}{{\mu }_s}\right)\left(\nabla ·{\overrightarrow{u}}_s\right)\overline{\overline{I}}\right]$$ (A2)

Granular bulk viscosity

$${\lambda }_s=\frac{4}{3}{\alpha }_s{\rho }_s{d}_s{g}_{0,ss}\left(1+e_{ss}\right){\left(\frac{{\Theta }_s}{\pi }\right)}^{1/2}$$ (A3)

Solid shear viscosity

$${\mu }_s={\mu }_{s,\mathrm{col}}+{\mu }_{s,\mathrm{kin}}+{\mu }_{s,\mathrm{fr}}$$ (A4)

Solid collision viscosity [32]

$${\mu }_{s,\mathrm{col}}=\frac{4}{3}{\alpha }_s^2{\rho }_s{d}_s{g}_{0,ss}\left(1+e_{ss}\right){\left(\frac{{\Theta }_s}{\pi }\right)}^{1/2}$$ (A5)

Solid kinetic viscosity

Syamlal-O`Brien [32] $${\mu }_{s,kin}=\frac{d_s{\rho }_s\sqrt{{\Theta }_s\pi }}{6\left(3-e_{ss}\right)}\left[1+\frac{2}{5}\left(1+e_{ss}\right)\left(3e_{ss}-1\right){\alpha }_s{g}_{0,ss}\right]$$ (A6)    Gidaspow [36] $${\mu }_{s,kin}=\frac{10{\rho }_s{d}_s\sqrt{{\Theta }_s\pi }}{96{\alpha }_s\left(1+e_{ss}\right)g_{0,ss}}{\left[1+\frac{4}{5}{g}_{0,ss}{\alpha }_s\left(1+e_{ss}\right)\right]}^2$$ (A7)

Solid frictional viscosity [27]

$${\mu }_s^{fr}=\frac{P_{s}sin\varnothing }{2\sqrt{I_{2D}}}$$ (A8) $$I_{2D}=\frac{1}{6}\left[{\left(D_{s,11}-D_{s,22}\right)}^2+{\left(D_{s,22}-D_{s,33}\right)}^2+{\left(D_{s,33}-D_{s,11}\right)}^2\right]$$ $$+D_{s,12}^2+D_{s,23}^2+D_{s,31}^2$$ (A9) $$D_{s,ij}=\frac{1}{2}\left(\frac{\partial {\overrightarrow{u}}_{s,i}}{\partial x_j}+\frac{\partial {\overrightarrow{u}}_{s,j}}{\partial x_i}\right)$$ (A10)

Solid phase pressure [26]

$$p_s={\alpha }_s{\rho }_s{\Theta }_s+2{\rho }_s\left(1+e_{ss}\right){\alpha }_s^2{g}_{0,ss}{\Theta }_s$$ (A11)

Radial distribution function [26,37]

$$g_{0,ss}=\frac{1}{1-{\alpha }_s}+\frac{3}{2}\frac{{\alpha }_s}{{(1-{\alpha }_s)}^2}+\frac{1}{2}\frac{{\alpha }_s^2}{{(1-{\alpha }_s)}^3}$$ (A12)

Table A1. Constitutive equations used in this study

Gas phase shear stress tensor

$${\overline{\overline{\tau }}}_g={\alpha }_g{\mu }_g\left[\nabla {\overrightarrow{u}}_g+{\left(\nabla {\overrightarrow{u}}_g\right)}^T-\frac{2}{3}\left(\nabla ·{\overrightarrow{u}}_g\right)\overline{\overline{I}}\right]$$ (A1)

Solid phase shear stress tensor

$${\overline{\overline{\tau }}}_s={\alpha }_s{\mu }_s\left[\nabla {\overrightarrow{u}}_s+{\left(\nabla {\overrightarrow{u}}_s\right)}^T+\left(\frac{2}{3}-\frac{{\lambda }_s}{{\mu }_s}\right)\left(\nabla ·{\overrightarrow{u}}_s\right)\overline{\overline{I}}\right]$$ (A2)

Granular bulk viscosity

$${\lambda }_s=\frac{4}{3}{\alpha }_s{\rho }_s{d}_s{g}_{0,ss}\left(1+e_{ss}\right){\left(\frac{{\Theta }_s}{\pi }\right)}^{1/2}$$ (A3)

Solid shear viscosity

$${\mu }_s={\mu }_{s,\mathrm{col}}+{\mu }_{s,\mathrm{kin}}+{\mu }_{s,\mathrm{fr}}$$ (A4)

Solid collision viscosity [32]

$${\mu }_{s,\mathrm{col}}=\frac{4}{3}{\alpha }_s^2{\rho }_s{d}_s{g}_{0,ss}\left(1+e_{ss}\right){\left(\frac{{\Theta }_s}{\pi }\right)}^{1/2}$$ (A5)

Solid kinetic viscosity

Syamlal-O`Brien [32] $${\mu }_{s,kin}=\frac{d_s{\rho }_s\sqrt{{\Theta }_s\pi }}{6\left(3-e_{ss}\right)}\left[1+\frac{2}{5}\left(1+e_{ss}\right)\left(3e_{ss}-1\right){\alpha }_s{g}_{0,ss}\right]$$ (A6)    Gidaspow [36] $${\mu }_{s,kin}=\frac{10{\rho }_s{d}_s\sqrt{{\Theta }_s\pi }}{96{\alpha }_s\left(1+e_{ss}\right)g_{0,ss}}{\left[1+\frac{4}{5}{g}_{0,ss}{\alpha }_s\left(1+e_{ss}\right)\right]}^2$$ (A7)

Solid frictional viscosity [27]

$${\mu }_s^{fr}=\frac{P_{s}sin\varnothing }{2\sqrt{I_{2D}}}$$ (A8) $$I_{2D}=\frac{1}{6}\left[{\left(D_{s,11}-D_{s,22}\right)}^2+{\left(D_{s,22}-D_{s,33}\right)}^2+{\left(D_{s,33}-D_{s,11}\right)}^2\right]$$ $$+D_{s,12}^2+D_{s,23}^2+D_{s,31}^2$$ (A9) $$D_{s,ij}=\frac{1}{2}\left(\frac{\partial {\overrightarrow{u}}_{s,i}}{\partial x_j}+\frac{\partial {\overrightarrow{u}}_{s,j}}{\partial x_i}\right)$$ (A10)

Solid phase pressure [26]

$$p_s={\alpha }_s{\rho }_s{\Theta }_s+2{\rho }_s\left(1+e_{ss}\right){\alpha }_s^2{g}_{0,ss}{\Theta }_s$$ (A11)

Radial distribution function [26,37]

$$g_{0,ss}=\frac{1}{1-{\alpha }_s}+\frac{3}{2}\frac{{\alpha }_s}{{(1-{\alpha }_s)}^2}+\frac{1}{2}\frac{{\alpha }_s^2}{{(1-{\alpha }_s)}^3}$$ (A12)

Table A1. Continued 

Diffussion coefficient of granular temperature

Syamlal-O`Brien [32]

$$k_{{\Theta }_s}=\frac{15d_s{\rho }_s{\alpha }_s\sqrt{{\Theta }_s\pi }}{4(41-33\eta )}$$ $$\left[1+\frac{12}{5}{\eta }^2(4\eta -3){\alpha }_s{g}_{0,ss}+\frac{16}{15\pi }(41-33\eta )\eta {\alpha }_s{g}_{0,s}\right]$$ (A13) where $$\eta =\frac{1}{2}\left(1+e_{ss}\right)$$ (A14)    Gidaspow [36] $$k_{{\Theta }_s}=\frac{150{\rho }_s{d}_s\sqrt{\left(\Theta \pi \right)}}{384\left(1+e_{ss}\right)g_{0,ss}}{\left[1+\frac{6}{5}{\alpha }_s{g}_{0,s}\left(1+e_s\right)\right]}^2$$ $$+2{\rho }_s{\alpha }_s^2{d}_s\left(1+e_{ss}\right)g_{0,ss}\sqrt{\frac{{\Theta }_s}{\pi }}$$ (A15)

Collision dissipation of energy [26]

$${\gamma }_{\Theta m}=\frac{12\left(1-e_{sS}^2\right)g_{0,ss}}{d_s\sqrt{\pi }}{\rho }_s{\alpha }_s^2{\Theta }_s^{3/2}$$ (A16)

Transfer of kinetic energy [36]

$${\phi }_{gs}=-3K_{gs}{\Theta }_s$$ (A17)

Table A1. Continued 

Diffussion coefficient of granular temperature

Syamlal-O`Brien [32]

$$k_{{\Theta }_s}=\frac{15d_s{\rho }_s{\alpha }_s\sqrt{{\Theta }_s\pi }}{4(41-33\eta )}$$ $$\left[1+\frac{12}{5}{\eta }^2(4\eta -3){\alpha }_s{g}_{0,ss}+\frac{16}{15\pi }(41-33\eta )\eta {\alpha }_s{g}_{0,s}\right]$$ (A13) where $$\eta =\frac{1}{2}\left(1+e_{ss}\right)$$ (A14)    Gidaspow [36] $$k_{{\Theta }_s}=\frac{150{\rho }_s{d}_s\sqrt{\left(\Theta \pi \right)}}{384\left(1+e_{ss}\right)g_{0,ss}}{\left[1+\frac{6}{5}{\alpha }_s{g}_{0,s}\left(1+e_s\right)\right]}^2$$ $$+2{\rho }_s{\alpha }_s^2{d}_s\left(1+e_{ss}\right)g_{0,ss}\sqrt{\frac{{\Theta }_s}{\pi }}$$ (A15)

Collision dissipation of energy [26]

$${\gamma }_{\Theta m}=\frac{12\left(1-e_{sS}^2\right)g_{0,ss}}{d_s\sqrt{\pi }}{\rho }_s{\alpha }_s^2{\Theta }_s^{3/2}$$ (A16)

Transfer of kinetic energy [36]

$${\phi }_{gs}=-3K_{gs}{\Theta }_s$$ (A17)