A Computational Study of the Influence of Drag Models and Heat Transfer Correlations on the Simulations of Reactive Polydisperse Flows in Bubbling Fluidized Beds

Abstract

In this work, the influence of gas–solid drag and heat transfer coefficient models on the prediction capacity of the Euler–Euler approach to simulate reactive bubbling fluidized bed flows is studied. Three different cases are considered, a non-reactive bidisperse bubbling fluidized bed flow (Case 1), and two reactive polydisperse flows in bubbling fluidized beds, one for biomass gasification (Case 2), and the other for biomass pyrolysis (Case 3). The Gidaspow, Syamlal–O’Brien, and BVK gas–solid drag models and the Gunn, Ranz–Marshall, and Li–Mason gas–solid heat transfer correlations are investigated. A Eulerian multiphase approach in a two-dimensional Cartesian domain is employed for the simulations. Computational results for the three cases are compared with experimental data from the literature. The results obtained here contribute to a better understanding of the impacts of such closure models on the prediction ability of the Euler–Euler approach to simulate reactive flows. The results indicate that, for the simulation of reactive flows in bubbling fluidized bed reactors, the kinetic modeling of the reactions has a global effect, which superposes with the influence of the drag and heat transfer coefficient models. Nevertheless, local parameters can be noticeably affected by the choice of the interface closure models. Finally, this work also identifies the models that lead to the best results for the cases analyzed here, and thus proposes the use of such selected models for gasification and pyrolysis processes occurring in bubbling fluidized bed reactors.

1. Introduction

The recent global awareness about the impact of the use of traditional fuels, such as petroleum and coal, has motivated a huge effort of the scientific and industrial communities all over the world to develop clean energy sources to supply the increasing demand for electric power. The pursuit for a near-zero emission source of energy has become a fundamental aspect to reach energetic sustainability, and to reduce the external dependence of petroleum in many countries [1]. Among all renewable resources, biomass is one of the most promising. One of the greatest advantages of biomass is that it contains carbon, as opposed to solar and wind sources, which permits a large range of applications, after an adequately thermochemical or biological treatment [2]. In fact, there are two possible ways to convert biomass in useful fuels: by thermochemical and by biochemical conversion routes [2,3]. The thermochemical route has several advantages in comparison with the biochemical route, such as the adequacy to produce diverse oxygenated and hydrocarbon fuels, reaction times that are several orders of magnitude shorter than biological processing, lower cost of catalysts, and the ability to recycle catalysts [3].

Pyrolysis and gasification are thermochemical conversion routes for biomass. The pyrolysis consists of a thermal degradation of the organic matter in the absence of oxidizing agents, yielding three groups of products: non-condensable gases, pyrolytic liquid (bio-oil/tar), and char [4]. In the gasification process, the thermal degradation occurs in the presence of a gasifying agent, which is usually air, steam, oxygen, or a mixture of these [2]. Both processes can be found in fixed and fluidized beds. When the processes occur in fluidized beds, inert or catalyzing solid particles are used to enhance mixing, and thus the heat and mass transfer inside the reactor. Further, there is a granular reactive flow occurring inside the reactor. As pyrolysis and gasification of biomass are significantly influenced by the heating rate of the fuel particles, the understanding of the fluid dynamics and heat transfer during both processes is essential to the development of these technologies. A very helpful and efficient tool to study the flow of such complex phenomena is the use of computational fluid dynamics (CFD).

In view of the recent advances in computer processing and numerical tools, many authors have used CFD to simulate the flow in fluidized beds in situations with and without chemical reactions. Biomass pyrolysis and gasification flows in fluidized beds have been simulated with CFD tools using both Euler–Euler [5,6,7,8,9,10,11,12] and Euler–Lagrange approaches [13,14,15,16,17,18]. Recently, Alobaid et al. [19] presented a very complete review concerning the current status of CFD modeling of fluidized beds for chemical and energy applications, which included the pyrolysis and gasification processes. The authors concluded that the Euler–Euler approach is indeed an efficient tool to study these processes.

Although those studies have successfully predicted the pyrolysis products leaving the reactor, in comparison to the experimental results, there is no quantitative estimate of the influence of closure models on the accuracy of the simulation results. The Euler–Euler approach demands a great number of closure models, which may have a significant impact on the predicted results. Among those closure models, the interface transport coefficients accounting for drag and heat transfer between the gas and solid phases are certainly of crucial importance. There are many works in the literature about drag models in bubbling fluidized [5,19,20,21,22,23,24] and circulating beds [10,19,25,26]. Unfortunately, these works are for monodisperse and non-reactive isothermal simulations only. In contrast, the flow inside a fluidized reactor during pyrolysis and gasification processes is polydisperse and reactive, with transport processes between the biomass and gas phases.

Generally speaking, the literature shows that, for a non-reactive bubbling fluidized bed, the Syamlal–O’Brien [27] and the Gidaspow [28] drag models may be considered appropriate [19], while for a non-reactive circulating fluidized bed, sub-grid models should be considered [19,25,26]. Brief comments on some recent studies of gas–solid closure models in Euler–Euler simulations of fluidized beds are outlined hereafter. Jia et al. [23] argued for polydisperse beds, a drag correlation originally based on bubble structure for monodisperse fluidized beds. The proposed model was compared with the BVK and Huilin–Gidaspow models, presenting better agreement with the experimental data for a bidisperse non-reactive fluidized bed. Varghese et al. [21] conducted 2D (two-dimensional) and 3D (three-dimensional) simulations of bubbling and turbulent fluidized beds, operating with monodisperse Geldart B particles using different gas–solid drag correlations. The results showed that the 2D and 3D simulations of bubbling fluidized beds using the Gidaspow model were able to adequately predict the mean void fraction, with an approximate error below 5% in comparison to the experimental data.

It should also be mentioned that the mass transfer between the biomass and gas phases may change the drag at their interface. However, simulation studies of reactive fluidized beds simply employ the Syamlal–O’Brien [5,6,7,9,10,11,12,29,30] drag model, with no further discussion. Only two works have been found that analyze the influence of the drag model in the simulation of a reactive bubbling bed. Xiong and Kong [5] tested the Syamlal–O’Brien, Gidaspow, and EMMS (Energy Minimum Multiscale Method) models for a biomass pyrolysis process. The results obtained by Xiong and Kong [5] suggested that the predicted pyrolysis products are significantly influenced by the drag models. Zhong et al. [24] included the effect of particle reduction in the Euler–Euler approach for biomass pyrolysis in a bubbling fluidized bed. The authors also compared axial results using three gas–solid drag models, Huilin–Gidaspow with Geldart B particles, EMMS with Geldart A particles, and a combination of the previous ones. The results between the three models were similar, except for the total mass of bio-mixture in the reactor, which was lower using the Huilin–Gidaspow model. However, the authors did not compare their computational results with the experimental data.

Now, little is known with respect to the influence of the heat transfer coefficient between the gas and solid phases on the prediction of pyrolysis products. Only Xiong and Kong [5] have studied this type of influence. The correlations provided by Gunn [31], Ranz–Marshall [32], and Li–Mason [33] were tested, and showed only a small impact of the model of interface heat transfer on the predicted products of wood pyrolysis. The correlation of Gunn is used in almost all works that were reviewed [10,11,12,24,34]. Only a few exceptions adopt the Ranz–Marshall correlation [7,35,36].

Based on the current literature, it is noted that there exists a need to better understand the influence of gas–solid interface coefficient models on reactive fluidized bed applications, such as biomass pyrolysis and gasification. The influence of kinetic models for these processes were evaluated in previous works [6,7]. In this work, the influence of drag and heat transfer correlations for gas–solids force and energy interactions are studied in three different cases: a non-reactive bidisperse bubbling fluidized bed (Case 1), and two reactive polydisperse flows in a bubbling fluidized bed, with one for biomass gasification (Case 2) and the other for biomass pyrolysis (Case 3). The Gidaspow, Syamlal–O’Brien, and BVK gas–solid drag models, and the Gunn, Ranz–Marshall, and Li–Mason gas–solid heat transfer correlations are investigated. Computational results are compared with the experimental data from the literature, and the most suitable models are identified for the situations studied here. The results presented here contribute to a better understanding of the impacts of such closure models on the prediction ability of the Euler–Euler approach to simulate reactive flows.

2. Mathematical Modeling

2.1. Governing Equations

The Euler–Euler approach [27] is employed to describe the multiphase reactive flow in the fluidized bed reactor, with the gas and all solid phases being assumed as interpenetrating continua. All phases are described using the Eulerian approach, and volume-averaged balance equations are written for mass, momentum, energy, and species conservation. A discussion about the assumptions adopted in the simulations of multiphase reactive flows are provided elsewhere [30,37].

The continuity equation for the gas phase is given by

$$\frac{\partial }{\partial t}\left({\alpha }_g{\rho }_g\right)+\nabla ·\left({\alpha }_g{\rho }_g{\mathbf{u}}_g\right)=\sum _{n=1}^{N_g}{R}_{gn}$$

(1)

where α_g, ρ_g, u_g, and N_g are the gas phase volume fraction, density, velocity, and species number, respectively, and R_gn is the net mass exchange between the solid phases and the nth species of the gas phase due to heterogeneous reactions.

The gas phase momentum equation can be written as [7]

$$\begin{gathered} \frac{\partial }{\partial t}\left({\alpha }_g{\rho }_g{\mathbf{u}}_g\right)+\nabla ·\left({\alpha }_g{\rho }_g{\mathbf{u}}_g{\mathbf{u}}_g\right)={\alpha }_g{\rho }_g\mathbf{g}-{\alpha }_g\nabla P_g+ \\ +\nabla ·\left({\mathsf{\tau }}_g+{\mathsf{\tau }}_g^t\right)+\sum _{m=1}^M{\beta }_{gm}\left({\mathbf{u}}_m-{\mathbf{u}}_g\right)+\sum _{m=1}^M{\mathsf{\psi }}_{gm} \end{gathered}$$

(2)

where P_g, ${\mathsf{\tau }}_g$, and ${\mathsf{\tau }}_g^t$ are the gas phase pressure, viscous stress tensor, and turbulent stress tensor, respectively, g is the acceleration of gravity, and ψ_gm accounts for the momentum generation due to the mass transfer induced by chemical reactions between the gas and mth solid phase. The κ-ε model detailed in Benyahia et al. [38] is employed to evaluate the turbulent stress tensor.

By neglecting variations in kinetic and potential energies, viscous dissipation, and volumetric expansion effects, the internal energy conservation equation is written as follows [7]:

$${\alpha }_g{\rho }_g{c}_{pg}\left(\frac{\partial T_g}{\partial t}+{\mathbf{u}}_g·\nabla T_g\right)=-\nabla ·\left({\mathbf{q}}_g+{\mathbf{q}}_g^t\right)+\sum _{m=1}^M{W}_{gm}-∆H_g$$

(3)

where T_g, c_pg, q_g, and ${\mathbf{q}}_g^t$ are the gas phase temperature, constant pressure specific heat, conductive heat flux, and turbulent heat flux, respectively, ΔH_g is the energy generation or consumption due to chemical reactions occurring in the gas phase, and W_gm is the heat transfer between the gas and mth solid phase.

The chemical species conservation equation for the gas phase is written as [7]

$$\frac{\partial }{\partial t}\left({\alpha }_g{\rho }_g{Y}_{gn}\right)+\nabla ·\left({\alpha }_g{\rho }_g{Y}_{gn}{\mathbf{u}}_g\right)=\nabla ·\left({\mathbf{J}}_{gn}+{\mathbf{J}}_{gn}^t\right)+R_{gn}$$

(4)

where Y_gn, J_gn, ${\mathbf{J}}_{gn}^t$, and R_gn represent the mass fraction, mass flux vector, turbulent mass flux vector, and production or consumption rate of the nth gas species, respectively.

The continuity equation for the mth solid phase, m = 1, …, N_m, is given by

$$\frac{\partial }{\partial t}\left({\alpha }_m{\rho }_m\right)+\nabla ·\left({\alpha }_m{\rho }_m{\mathbf{u}}_m\right)=\sum _{n=1}^{N_m}{R}_{mn}$$

(5)

where α_m, u_m, and ρ_m are the mth solid phase volume fraction, velocity, and density, respectively. The term R_mn is the net mass exchange rate between the mth solid phase and gas phase due to heterogeneous reactions.

The momentum equation for the mth solid phase, m = 1, …, N_m, is given by [7]

$$\begin{gathered} \frac{\partial }{\partial t}\left({\alpha }_m{\rho }_m{\mathbf{u}}_m\right)+\nabla ·\left({\alpha }_m{\rho }_m{\mathbf{u}}_m{\mathbf{u}}_m\right)=-{\alpha }_m\nabla P_g+\nabla ·{\mathsf{\tau }}_m+{\alpha }_m{\rho }_m\mathbf{g}+{\beta }_{gm}\left({\mathbf{u}}_m-{\mathbf{u}}_g\right) \\ -\sum _{l=1}^M{\beta }_{lm}\left({\mathbf{u}}_m-{\mathbf{u}}_l\right)+{\mathsf{\psi }}_{gm}+\sum _{l=1}^M{\mathsf{\psi }}_{lm} \end{gathered}$$

(6)

where τ_m is the stress tensor of the mth solid phase; β_ml and ψ_lm are, respectively, the drag coefficient and mass transfer induced source of momentum between the solid phases l and m.

In order to evaluate the solid stress tensor, the kinetic theory of granular flow (KTGF) is considered. KTGF assumes the solid stress tensor to be a function of the granular temperature, Θ_m, which is obtained using an algebraic equation, given as [27]

$${\Theta }_m=\left\{\frac{K_{1m}{\alpha }_{m}tr\left({\mathbf{D}}_m\right)+\sqrt{K_{1m}^2{tr}^2\left({\mathbf{D}}_m\right){\alpha }_m^2+4K_{4m}{\alpha }_m\left[K_{2m}{tr}^2\left({\mathbf{D}}_m\right)+2K_{3m}tr\left({\mathbf{D}}_m^2\right)\right]}}{2{\alpha }_m{K}_{4m}}\right\}$$

(7)

where D_m is the rate of strain tensor, and the expressions for K_1m, K_2m, K_3m, and K_4m are g in the literature [27].

By neglecting the heat transfer between different solid phases, the energy conservation equation for the solid phase m can be written as [6,7,27]

$${\alpha }_m{\rho }_m{c}_{pm}\left(\frac{\partial T_m}{\partial t}+{\mathbf{u}}_m·\nabla T_m\right)=-\nabla ·{\mathbf{q}}_m-W_{gm}-∆H_m$$

(8)

where T_m, c_pm, and q_m are, respectively, the temperature, constant pressure specific heat, and conductive heat flux for the mth solid phase, and ΔH_m is the heat generated/consumed due to the chemical reactions occurring in the solid phase m. It is assumed that the solid phase obeys the Fourier law [27].

The chemical species conservation equation for the solid phase m is written as [7]

$$\frac{\partial }{\partial t}\left({\alpha }_m{\rho }_m{Y}_{mn}\right)+\nabla ·\left({\alpha }_m{\rho }_m{Y}_{mn}{\mathbf{u}}_m\right)=\nabla ·{\mathbf{J}}_{mn}+R_{mn}$$

(9)

where Y_mn, J_mn, and R_mn represent the mass fraction, mass flux vector, and production or consumption rate of the nth species in the mth solid phase, respectively.

It must also be noted that the Soret effect is not considered here, and the diffusivity coefficient of gaseous species is evaluated using the Bird model. Furthermore, the source terms related to the chemical reactions are evaluated according to the appropriate kinetic modeling.

2.2. Interface Transport Coefficients Correlations

The momentum and energy conservation equations need closure models to describe the interaction between the phases, and this is performed with the interface transport coefficients. The transfer of momentum between phases includes different interaction effects. Among the most important are the (i) drag forces due to the differences in phase velocities, (ii) buoyancy forces due to the pressure gradient in the gas phase, and (iii) the momentum transfer due to the mass transfer. The buoyancy forces are considered in Equation (2), while the momentum transfer due to mass transfer is associated with the assumed heterogeneous reactions. The drag forces are modeled using different equations to evaluate the drag coefficient ${\beta }_{gm}$.

Here, the influence of different drag coefficient correlations on granular flow prediction for reactive and non-reactive systems is studied by adopting three different drag models from the literature. Two of them were formulated for monodisperse flows (Gidaspow [28] and Syamlal–O’Brien [27]), and the other one was formulated for polydisperse flows (BVK [39]). The three correlations are described in the following paragraphs.

• Syamlal–O’Brien: It is based on the terminal velocity of a fluidized bed. It was derived for a single spherical particle, and modified for inclusion in Euler–Euler models for granular flows. The gas–solid drag coefficient is given by

  $${\beta }_{gm}=\frac{3}{4}{C}_{Dm}\frac{{\epsilon }_m{\epsilon }_g{\rho }_g}{u_{rm}^2{d}_{pm}}\left|{\mathbf{u}}_g-{\mathbf{u}}_m\right|$$

  (10a)

  where $d_{pm}$ is the mean diameter of the particles of the solid phase m. $C_D$ and u_rm are, respectively, the drag coefficient of one particle and the terminal velocity, which are evaluated as

  $$C_{Dm}={\left(0.063+4.8\sqrt{u_{rm}/{\mathrm{Re}}_m}\right)}^2$$

  (10b)

  $$u_{rm}=0.5\left[\sqrt{{\left(0.06{\mathrm{Re}}_m\right)}^2+0.12{\mathrm{Re}}_m\left(2b-a\right)+a^2}+a-0.06{\mathrm{Re}}_m\right]$$

  (10c)

  Re_m is the Reynolds number. The coefficients a and b are given, respectively, as

  $$a={\epsilon }_g^{4.14}$$

  (10d)

  $$b=\left\{\begin{array}{c}0.8{\epsilon }_g^{1.28},{\epsilon }_g\le 0.85\\ {\epsilon }_g^{2.65},{\epsilon }_g>0.85\end{array}\right.$$

  (10e)
• Gidaspow: It is a combination of the Ergun equation (ε_g ≥ 0.8)—obtained for a dense fixed bed and based in the pressure drop in the bed—and Wen-Yu equation (ε_g < 0.8)—which was developed for a fluidized bed at minimum fluidization condition [28]. In this model, the gas–solid drag coefficient is then given by

  $${\beta }_{gm}=\left\{\begin{array}{c}150\frac{{\epsilon }_m^2{\mu }_g}{d_{pm}^2{\epsilon }_g}+1.75\frac{{\epsilon }_m{\rho }_g}{d_{pm}}\left|{\mathbf{u}}_g-{\mathbf{u}}_m\right|,{\epsilon }_g<0.8\\ \frac{3}{4}{C}_{Dm}\frac{{\epsilon }_m{\rho }_g}{d_{pm}}\left|{\mathbf{u}}_g-{\mathbf{u}}_m\right|{\epsilon }_g^{-2.65},{\epsilon }_g\ge 0.8\end{array}\right.$$

  (11a)

  where the drag coefficient of one particle, C_D, is written as

  $$C_{Dm}=\left\{\begin{array}{c}\frac{24}{{\epsilon }_g{\mathrm{Re}}_m}\left[1+0.15{\left({\epsilon }_g{\mathrm{Re}}_m\right)}^{0.687}\right],{\mathrm{Re}}_m<1000\\ 0.44,{\mathrm{Re}}_m\ge 1000\end{array}\right.$$

  (11b)
• BVK: It was developed by Beetstra et al. [39] from Lattice–Boltzmann simulations for polydisperse flows. The gas–solid drag coefficient in this model is given as

  $${\beta }_{gm}=180\frac{{\mu }_g}{d_{pm}^2}\frac{{\left(1-{\epsilon }_g\right)}^2}{{\epsilon }_g}+18\frac{{\mu }_g}{d_{pm}^2}{{\epsilon }_g}^3\left(1-{\epsilon }_g\right)\left(1+1.5\sqrt{1-{\epsilon }_g}\right)$$

  (12)

  The heat transfer between gas and solid phases, W_gm, is expressed by

  $$W_{gm}={\gamma }_{gm}\left(T_m-T_g\right)$$

  (13)

  The term ${\gamma }_{gm}$ represents the heat transfer coefficient between phases, and is given by the following formula:

  $${\gamma }_{gm}=\frac{6k_g{\epsilon }_m{\mathrm{Nu}}_m}{d_{pm}^2}$$

  (14)

  where k_g is the gas phase thermal conductivity. To obtain the Nusselt number of the heat transfer between phases, three different correlations were studied: Gunn [31], Ranz–Marshall [32], and Li–Mason [33].
• Gunn: It includes the effect of the gas volumetric fraction in its equations. According to Di Natale et al. [40], the Gunn model is a robust expression, and allows for an accurate evaluation of the gas–solid heat transfer coefficient. In the Gunn model, the Nusselt number is computed as

  $${\mathrm{Nu}}_m=\left(7-10{\epsilon }_g+5{\epsilon }_g^2\right)\left(1+0.7{\mathrm{Re}}_m^{0.2}{\mathrm{Pr}}^{1/3}\right)+\left(1.33-2.4{\epsilon }_g+1.2{\epsilon }_g^2\right){\mathrm{Re}}_m^{0.7}{\mathrm{Pr}}^{1/3}$$

  (15)
• Ranz–Marshall: It was derived to describe the heat exchanged between a single particle and a homogeneous fluid. Nevertheless, it is shown that the model can be used to evaluate the Nusselt number in fluidized beds with Re_m > 100 [40,41]. In this model, the Nusselt number is given as

  $${\mathrm{Nu}}_m=2+0.6{\mathrm{Re}}_m^{1/2}{\mathrm{Pr}}^{1/3}$$

  (16)
• Li–Mason: It includes the void fraction, and assumes three different expressions for the Nusselt number in terms of the Reynolds number. Its mathematical expressions are given as

  $${\mathrm{Nu}}_m=\left\{\begin{array}{c}2+0.6{\epsilon }_g^{3.5}{\mathrm{Re}}_m^{1/2}{\mathrm{Pr}}^{1/3},{\mathrm{Re}}_m\le 200\\ 2+0.5{\epsilon }_g^{3.5}{\mathrm{Re}}_m^{1/2}{\mathrm{Pr}}^{1/3}+0.02{\epsilon }_g^{3.5}{\mathrm{Re}}_m^{0.8}{\mathrm{Pr}}^{1/3},200<{\mathrm{Re}}_m\le 1500\\ 2+0.000045{\epsilon }_g^{3.5}{\mathrm{Re}}_m^{0.8},{\mathrm{Re}}_m>1500\end{array}\right.$$

  (17)

  For all equations above, the Reynolds and Prandtl numbers are given, respectively, by

  $${\mathrm{Re}}_m=\frac{d_{pm}\left|{\mathbf{u}}_g-{\mathbf{u}}_m\right|{\rho }_g}{{\mu }_g}$$

  (18a)

  $$\mathrm{Pr}=\frac{c_{pg}{\mu }_g}{k_g}$$

  (18b)

3. Studied Cases

In order to study the influence of the gas–solid drag and heat transfer coefficients on the prediction of granular flows in reactive bubbling fluidized beds, three cases are considered here, as follows: a non-reactive bidisperse flow (Case 1), and two reactive flows, one with wood gasification (Case 2) and one with sugarcane bagasse pyrolysis (Case 3). The goal is to evaluate the effects of the models on the flow behavior in different situations.

For the cases analyzed here, a 2D computational domain is used to describe the bubbling fluidized bed reactor. In the gasification and pyrolysis cases, the literature shows that such 2D approach is valid if the superficial velocity of the fluidizing agent and the ratio of biomass to gas mass flow rates are the same as for the 3D real domain [6,7,29]. The computational modeling was completed using the MFIX (Multiphase Flow with Interphase eXchanges) code [42]. A modified SIMPLE algorithm [43] is adopted to solve the governing equations. Furthermore, in the three studied cases, a Total Variation Diminishing (TVD) scheme is used for the spatial discretization of the convective terms, using the Smart [44] flux limiter. In order to reduce computational costs, an adaptive time step is used. The marching scheme for the time integration of the governing equations is the implicit Euler scheme. Residue values of 10⁻³ are adopted for continuity and momentum equations, and values of 10⁻⁴ are used for species and energy equations.

The Euler–Euler closure models employed in the three cases investigated here are listed in Table 1. The shear stress model assumes that, for the volume fraction of the mth phase smaller than critical packing, ${\epsilon }_m^{*}$, the granular plastic model of Schaeffer is adopted, while for the volume fraction of the mth phase higher than or equal to the critical packing, the KTGF formulation from Agrawal is used. Further details of the three cases investigated here are given in Section 3.1, Section 3.2 and Section 3.3.

3.1. Case 1: Non-Reactive Bidisperse Bubbling Fluidized Bed

Due to the lack of solids’ distribution data in reactive bubbling fluidized bed flows, a non-reactive flow is first analyzed. The flow described in Zhong et al. [45] is chosen because the authors presented experimental data for the solids’ volume fraction along the bed. In [45], a computational study was developed to investigate the influence of specularity, particle–particle restitution, and particle–wall coefficients on the predicted results of the CFD model based on the Euler–Euler approach.

Table 1. Closure models adopted in the cases studied in this work.

	Case 1	Case 2	Case 3
Flow	Non-reactive (bidisperse)	Reactive (wood gasification)	Reactive (sugarcane bagasse pyrolysis)
$\mathrm{Critical}\mathrm{packing}({\epsilon }_m^{*}$)	Fedors–Landel correlation	0.35	0.42
Shear stress model	${\epsilon }_m>{\epsilon }_m^{*}$ → Schaeffer model ${\epsilon }_m\le {\epsilon }_m^{*}$ → Agrawal model
Granular energy	Algebraic equation, Equation (7)
Fml	Syamlal Model [46]
Fgm	Syamlal–O’Brien, Gidaspow, and BVK
Wgm	--	Gunn, Ranz–Marshall, and Li–Mason
km	--	Constant and equal to the particle of mth phase
kg	--	Bird equation	Non-condensable gas: 2.577 × 10−2 W/m K Tar: 5.63 × 10−2 W/m K
μg	1.83 × 10−5 Pa·s	Sutherland equation	3 × 10−5 Pa s

Table 1. Closure models adopted in the cases studied in this work.

	Case 1	Case 2	Case 3
Flow	Non-reactive (bidisperse)	Reactive (wood gasification)	Reactive (sugarcane bagasse pyrolysis)
$\mathrm{Critical}\mathrm{packing}({\epsilon }_m^{*}$)	Fedors–Landel correlation	0.35	0.42
Shear stress model	${\epsilon }_m>{\epsilon }_m^{*}$ → Schaeffer model ${\epsilon }_m\le {\epsilon }_m^{*}$ → Agrawal model
Granular energy	Algebraic equation, Equation (7)
Fml	Syamlal Model [46]
Fgm	Syamlal–O’Brien, Gidaspow, and BVK
Wgm	--	Gunn, Ranz–Marshall, and Li–Mason
km	--	Constant and equal to the particle of mth phase
kg	--	Bird equation	Non-condensable gas: 2.577 × 10−2 W/m K Tar: 5.63 × 10−2 W/m K
μg	1.83 × 10−5 Pa·s	Sutherland equation	3 × 10−5 Pa s

The situation analyzed here is a bidisperse bubbling fluidized flow constituted by particles of glass and polystyrene. The particles are fluidized by a current of air that enters through the bottom. The experimental setup which was considered was a cylindrical column with 184 mm in diameter and 400 mm in height. The 2D domain adopted to simulate the flow is illustrated in Figure 1. The properties of the solid particles and air, and the assumed initial and boundary conditions are presented, respectively, in

Table 2. Properties of the particles for Case 1.

	ρk [kg/m3]	dpm [mm]	umf [cm/s]	μg [Pa·s]	${\mathit{\epsilon }}_{\mathit{m}}^{*}$	emm	Φ	Cf
Glass	2476	0.116	1.80	----	0.564	0.99	55°	0.1
Polystyrene	1064	0.275	4.00	----	0.574	0.99	55°	0.1
Gas	1.0	----	----	1.83 × 10−5	----	----	----	----

and

Table 3. Initial and boundary conditions.

		Pg	ug	um	${\mathit{\epsilon }}_{\mathit{g}}$	${\mathit{\epsilon }}_{\mathit{m}}$
Initial Conditions	Bed	1 atm	5 cm/s	0 cm/s	0.41	0.295
	Freeboard	1 atm	7 cm/s	0 cm/s	1	0
Boundary Conditions	Gas inlet	1 atm	3.84 cm/s	0 cm/s	1	0
	Wall	1 atm	No slip	No slip	1	0

. In Table 2, e_mm is the particle–particle restitution coefficient, ${\epsilon }_m^{*}$ is the critical packing, Φ is the internal friction angle, and C_f is the friction coefficient between phases. The mesh constructed has the same dimensions of the one used in the reference work [45], which consists of cells equally spaced by 5 mm in the x and y directions.

The closure models adopted in this work for Case 1 are described in Table 1. As one may note by viewing the table, the critical packing factor for Case 1 is evaluated through the Fedors–Landel model for bidisperse flows, which is given as

$${\epsilon }_s^{*}=\left\{\begin{array}{c}\begin{array}{c}\left[{\epsilon }_m^{*}-{\epsilon }_l^{*}+\left(1-\sqrt{\frac{d_{pl}}{d_{pm}}}\right)\left(1-{\epsilon }_m^{*}\right){\epsilon }_l^{*}\right]\\ \times \left[{\epsilon }_m^{*}+\left(1-{\epsilon }_m^{*}\right){\epsilon }_l^{*}\right]\frac{Y_m}{{\epsilon }_m^{*}}+{\epsilon }_l^{*}\end{array},\mathrm{if}{Y}_m\le \frac{{\epsilon }_m^{*}}{{\epsilon }_m^{*}+\left(1-{\epsilon }_m^{*}\right){\epsilon }_l^{*}}\\ \begin{array}{c}\left(1-\sqrt{\frac{d_{pl}}{d_{pm}}}\right)\left[{\epsilon }_m^{*}+\left(1-{\epsilon }_m^{*}\right){\epsilon }_l^{*}\right]\\ \times \left(1-Y_m\right)+{\epsilon }_m^{*}\end{array},\mathrm{if}{Y}_m>\frac{{\epsilon }_m^{*}}{{\epsilon }_m^{*}+\left(1-{\epsilon }_m^{*}\right){\epsilon }_l^{*}}\end{array}\right.$$

(19a)

where ${\epsilon }_m^{*}$ and ${\epsilon }_l^{*}$ are the critical packing factor of phases m and l, respectively, and Y_m is the mass fraction of solid phase m, which is defined as

$$Y_m=\frac{{\epsilon }_m^{*}}{{\epsilon }_m^{*}+{\epsilon }_l^{*}}$$

(19b)

3.2. Case 2: Biomass Gasification in a Bubbling Fluidized Bed

The gasifier studied in [29] is a lab scale reactor, in which the bed has a height of 600 mm and an internal diameter of 95 mm, while its freeboard is 450 mm high and has an internal diameter of 134 mm. A conical region of 50 mm connects the bed and the freeboard regions. There is a distribution plate located at the bottom of the reactor.

The gasifier domain is illustrated in Figure 2. In the considered domain, air enters through the bottom, and wood is inserted through a lateral port with 50 mm of diameter. The yielded gases leave the reactor at the top.

A noticeable difference between the biomass gasification process assumed in the reference work [29] and other works in the literature is the use of a bed constituted only by reactive particles, e.g., wood and char, without any inert material. In order to consider the reduction of char particles, three solid phases are assumed, i.e., wood, char 1 and char 2, where char 1 and char 2 have different diameters.

The governing equations described in Section 2 and the closure models presented in Table 1 are adopted for Case 2. The modeling of drying, pyrolysis, oxidation, and reduction reactions is made following the reference work [29], and the chemical reactions and respective kinetic parameters adopted are detailed in Appendix A. The initial and boundary conditions are also chosen following [29], and the values assumed for each quantity are summarized in

Table 4. Boundary conditions assumed for Case 2.

Wood Inlet		Air Inlet		Walls
uwood	0.035 cm/s	uair	25 cm/s	Tbed	970 K
Twood	423.15 K	Tair	670 K	Tfreeb	570 K
uH2O	8.406 cm/s	Pg	1 atm	ug	No slip
TH2O	423.15 K			um	No slip
YH2O	1
Pg	1 atm

and

Table 5. Initial conditions assumed for Case 2.

Initial Conditions for Case 2
Tg	1020 K	Tchar1	1020 K	Tchar2	1020 K
Pg	1 atm	uchar1	0.0 cm/s	Vchar2	0.0 cm/s
Vg	0.0 cm/s	εchar1	0.325	εchar2	0.325
YN2	1	εwood	0.0
εg	0.35	hbed	35 cm

. Prescribed temperatures are assumed at the gasifier walls. In the bed region, a prescribed temperature of 970 K is adopted, while in the freeboard region, a value of 570 K is assumed. These values were suggested by Gerber et al. [29] after a series of tests with different prescribed temperatures at the walls. Furthermore, it should be mentioned that, as the drying step was assumed instantaneous, it is established that 10% of the wood entering the gasifier is formed by steam.

As shown in Table 5, the reactor is initially filled with nitrogen (N₂) at rest and at atmospheric pressure, and contains a bed with height of 35 cm. The properties of the solid and gas phases adopted here are presented in

Table 6. Properties and physical parameters assumed in Case 2.

	ρk [kg/m3]	dpm [mm]	cp [J/kg K]	k [W/m K]	μ [Pa·s]	emm	Φ	Cf
Wood	585	4.0	2380	0.158	See Table 1	0.8	30°	0.1
Char 1	450	2.0	1600	0.107	See Table 1	0.8	30°	0.1
Char 2	450	1.0	1600	0.107	See Table 1	0.8	30°	0.1
Gas	Ideal gas	----	Ideal gas mixture with NASA polynomial	Bird equation	Sutherland equation	----	----	----

. The properties of the solid phases follow the ones from the reference work [29].

A convergence grid was constructed for five different meshes: 60 × 160, 20 × 240, 60 × 240, 80 × 320 and 35 × 196. The mesh with 35 × 196 was not equally spaced, with mesh elements larger in the bed region and smaller in the freeboard region. This mesh generated converged results, as well as the more refined one, but with a smaller computational cost. The details of the convergence grid are not given here, and can be found elsewhere [37].

3.3. Case 3: Biomass Pyrolysis in a Bubbling Fluidized Bed

Case 3 considers the geometry given in [47] for the sugarcane bagasse pyrolysis. The reactor has a height of 64 cm with a diameter that varies with the height, as seen in Figure 3. The biomass enters the reactor with a temperature of 50 °C and at a mass flow rate of 0.9 kg/h through a lateral port with 52 mm in diameter, and located 168.5 mm from the bottom. Nitrogen is used as the fluidizing agent, which enters the reactor through the bottom with a temperature of 773 K and a velocity of 37.19 cm/s.

The governing equations described in Section 2 and the closure models presented in Table 1 are adopted for Case 3. Drying is assumed instantaneous, and the sugarcane bagasse pyrolysis is modeled using the RAN kinetic scheme [6,48]. Further details about the kinetic model of pyrolysis chemical reactions are given in Appendix B. The initial and boundary conditions are summarized in

Table 7. Boundary conditions assumed for Case 3.

Boundary Conditions for Case 3
Sugarcane Bagasse Inlet			Nitrogen Inlet		Walls
$\dot{m}$biomass	0.0568 g/s		uN2	37.19 cm/s	Twalls		779 K
Tbiomass	323 K		TN2	773 K	ug		No slip
$\dot{m}$N2	0.0265 cm/s		Pg	1 atm	um		No slip
TN2	773 K
YN2	1
Pg	1 atm
Initial Conditions for Case 3
	Bed				Freeboard
	Gas	Sugarcane Bagasse	Sand	Char	Gas	Sugarcane Bagasse	Sand	Char
εm	0.4119	0.0001	0.58	0	0.999	0.001	0	0
|um| [cm/s]	0	0	0	0	0	0	0	0
Tm [K]	773 K				773 K
YN2	1				1

. A prescribed temperature of 773 K is assumed at the reactor’s walls. This value is based on a previous study [49]. In the experimental work of Hugo [47], the bed height is not informed; however, Xiong et al. [35] reported that the bed height has no effects on the distribution of pyrolysis products. Therefore, a bed height of 22 cm is considered. The no slip condition is enforced for both phases in this work. The properties of solid particles and the physical parameters assumed for Case 3 are shown in

Table 8. Properties and physical parameters assumed for Case 3.

	ρk [kg/m3]	dpm [mm]	cp [J/kg K]	k [W/m K]		μ [Pa·s]	emm	Φ	Cf
Sugarcane bagasse	200	2.0	1760	0.1		See Table 1	0.97	55°	0.1
Char	200	0.012	1100	0.1		See Table 1	0.97	55°	0.1
Sand	2650	0.5	800	0.27		See Table 1	0.97	55°	0.1
Gas	Ideal gas	----	Ideal gas mixture with NASA polynomial	NCG	2.577 × 10−2	3 × 10−5	----	----	----
				Tar	5.63 × 10−2

, where the values for specific heats, dynamic viscosity, and thermal conductivity are taken from [50].

A grid convergence test with three different mesh sizes was realized, showing that the mesh size of 38 × 128 provides a good balance between accuracy and computational time [6]. Therefore, this mesh size is employed here.

4. Results and Discussion

Results for the three cases described in Section 3.1, Section 3.2 and Section 3.3 are presented and discussed in Section 4.1, Section 4.2, and Section 4.3, respectively. All the results presented here are temporal average values for the last 20 s of the simulations in statistically stationary regimes.

4.1. Case 1

Results for the non-reacting bidisperse bubbling fluidized bed (Case 1) are discussed first. The mass fractions for the solid phases obtained in the simulations, when the interface gas–solid drag coefficient is varied, are compared to the experimental data from the reference work of Zhong et al. [45]. Figure 4 shows the axial profiles of the glass solid phase mass fraction obtained using the gas–solid drag coefficient from the Gidaspow, BVK, and Syamlal–O’Brien models. The experimental data reported by Zhong et al. [45] are also presented in Figure 4. For the conditions assumed for Case 1, the reference work [45] reports a segregation in the bed, with glass particles concentrated at the bottom of the bed. Figure 4 shows that Gidaspow and Syamlal–O’Brien’s correlations are able to adequately predict the segregation behavior of the bed. However, the BVK correlation is not able to describe the segregation behavior since it results in a more homogeneous distribution along the bed.

Figure 5 shows the radial profiles of glass phase mass fraction at different bed heights. The radial profiles obtained with the Gidaspow and Syamlal–O’Brien correlations present a behavior similar to the experimental one, with a higher concentration of glass particles at the center of the bed, and a small concentration near the walls. The Gidaspow and Syamlal–O’Brien correlations are also able to predict the symmetric profiles of glass particle distributions shown in the experimental data. Only the BVK correlation results in a non-symmetrical profile for dimensionless heights of 0.07 and 0.21. Nevertheless, despite the different radial profile, Figure 5 shows that the glass phase mass fraction predicted using the BVK gas–solid drag correlation is satisfactory compared to the experimental data in [45]. Based on the analysis of Figure 4 and Figure 5, it is apparent that the Gidaspow gas–solid drag correlation leads to the best prediction of the glass phase mass fraction. As a matter of fact, the results obtained with the Gidaspow correlation display the same trends in the axial and radial profiles observed in the experimental data [45], and the values of glass phase mass fraction are closer to the experimental data than those of the Syamlal–O’Brien gas–solid drag correlation.

Figure 6 and Figure 7 present the radial profiles of the glass and polystyrene velocities, respectively, obtained with the gas–solid drag correlations investigated here. It is known from the literature [45] that, in a bubbling fluidized bed, the solid particles move upwards at the center of the bed, and downwards near the walls, which implies annular radial profiles along the bed with positive axial velocity at the center, and negative axial velocity near the walls. The expected behavior is captured when the gas–solid drag correlations of Gidaspow and Syamlal–O’Brien are used at different bed heights, but not when the BVK correlation is used. Furthermore, Figure 6 and Figure 7 also reveal similar radial profiles of solid velocities for Gidaspow and Syamlal–O’Brien correlations, with the former leading to the prediction of higher velocities at the center of the bed and near the walls.

Results presented here for Case 1 indicate that different choices for the gas–solid drag correlation can cause significant changes in the flow behavior of non-reactive bubbling fluidized beds with different kinds of solid particles. For Case 1, the Gidaspow and Syamlal–O’Brien gas–solid drag correlations lead to similar behavior, with the Gidaspow correlation yielding the glass phase mass fraction closer to the experimental value. It is fair to state that the Gidaspow gas–solid drag correlation leads to the most satisfactory results for the situation studied in Case 1.

4.2. Case 2

Here, the gas–solid drag and heat transfer correlations are analyzed for wood gasification in a bubbling fluidized bed, which is described in Section 3.2. The experimental data from the literature [29] for yielded gas and gas temperature along the reactor are used to evaluate the results obtained with the different correlations studied. The prediction capacities of the simulations are compared by calculating the root mean square (RMS) to evaluate the yielded gas composition, and to obtain the relative difference (Dif_R) to evaluate the temperature of the gases leaving the reactor. In

Table 9. Experimental ranges of molar fractions and gas temperature at the reactor’s outlet reported in [29], and the corresponding mean values.

	Minimum Value	Maximum Value	Mean Value
XCO [%]	13.0	21.0	17.0
XCO2 [%]	13.0	17.0	15.0
XH2 [%]	7.0	11.0	9.0
XCH4 [%]	2.0	6.0	4.0
XN2 [%]	48.0	52.0	50.0
Tout [°C]	759.0	921.0	840.0

,

Table 10. Molar fractions of gas species and gas temperature at the reactor’s outlet for different gas–solid drag correlations.

Correlation	XN2 [%]	XH2 [%]	XCO [%]	XCH4 [%]	XH2O [%]	XCO2 [%]	RMS [%] a	Tout [K]	DifR [%] b
Gidaspow	50.4	9.2	17.9	3.0	3.8	13.5	0.9	906.4	7.9
Syamlal–O’Brien	49.2	10.0	21.5	2.8	2.6	11.6	2.6	871.7	3.8
BVK	50.9	9.1	17.2	2.9	3.8	13.9	0.8	887.5	5.7

^a Root mean square for the differences of gas molar fractions obtained in the simulations and the mean values of experimental ranges; ^b relative difference between the gas outlet temperature provided by the simulations and the mean value of the experimental range.

and

Table 11. Molar fractions of gas species and gas temperature at the reactor’s outlet for different gas–solid heat transfer correlations.

Correlation	XN2 [%]	XH2 [%]	XCO [%]	XCH4 [%]	XH2O [%]	XCO2 [%]	RMS [%] a	Tout [K]	DifR [%] b
Gunn	50.4	9.2	17.9	3.0	3.8	13.5	0.9	906.4	7.9
Ranz–Marshall	51.0	9.4	16.9	2.8	3.5	14.1	0.8	854.0	1.7
Li–Mason	50.7	9.1	16.6	3.0	4.1	14.2	0.7	901.6	7.3

^a Root mean square for the differences of gas molar fractions obtained in the simulations and the mean values of experimental ranges; ^b relative difference between the gas outlet temperature given by the simulations and the mean value of the experimental range.

, presented in this subsection, both the RMS and the Dif_R are calculated in terms of the center of the experimental range available in the reference work of Gerber et al. [29]. The experimental ranges presented in [29], and the corresponding mean values to calculate the RMS and Dif_R, are shown in Table 9.

First, the influences of gas–solid drag correlations are analyzed. Table 10 presents the results for molar fractions and gas temperatures leaving the reactor obtained with the Gidaspow, Syamlal–O’Brien, and BVK correlations. Calculated RMS and Dif_R are also shown. One may note from viewing Table 10 that the accuracies in the predicted yielded gas with the Gidaspow and BVK gas–solid correlations are very close, with RMS values of 0.9% and 0.8%, respectively. The gas–solid correlation of Syamlal–O’Brien leads to a slightly higher RMS value, which is equal to 2.6%. Concerning the gas temperatures leaving the reactor, the smallest and largest relative differences are found using the Syamlal–O’Brien (Dif_R = 3.8%) and Gidaspow (Dif_R = 7.9%) correlations, respectively.

Figure 8 shows the axial profiles of the gas temperatures provided by the simulations using the three gas–solid drag correlations studied here, as well as the experimental data reported in [29]. The experimental data in Figure 8 verifies that the gas temperature along the particle bed is almost constant, and that it noticeably drops at the end of the bed. This behavior is correctly predicted in the simulations with the three gas–solid drag models studied here, although the predicted height where the temperature drop occurs is smaller than the experimental value. The analysis illustrated in Figure 8 also shows that the Syamlal–O’Brien correlation results in a gas temperature axial profile which is closer to the experimental data, in comparison to those obtained with the Gidaspow and BVK correlations. One may also note from viewing Figure 8 that the three gas–solid drag correlations predict a sudden increase in gas temperature at the reactor’s bottom, followed by a temperature drop. On the one hand, the temperature rise is caused by the thermal energy released by the combustion reactions, which occur as soon as the air enters the gasifier. On the other hand, the temperature drop is caused by the wood entering the reactor at a low temperature, reducing the gas temperature near the wood-feeding position.

In order to evaluate the solids’ distribution along the reactor, the axial profiles for the volumetric fraction of all phases obtained with the three gas–solid drag correlations are analyzed. The axial profiles for the volumetric fraction of the gas phase are presented in Figure 9. It is clear that the Syamlal–O’Brien correlation results in higher concentrations of particles at heights below 42 cm than those predicted using the Gidaspow and BVK correlations. As a consequence, in the simulation using the Syamlal–O’Brien correlation, the height of the particle bed is smaller.

Figure 10 presents the axial profiles for the volumetric fractions of solid phases for char 1 (Figure 10a) and char 2 (Figure 10b). In Figure 10, the Syamlal–O’Brien correlation is used for the simulation, and a higher concentration of char at the bed region until a height of 42 cm is obtained. The results presented in Figure 9 and Figure 10 suggest that, for Case 2, the three drag correlations studied here lead to a similar behavior of the volumetric fractions of the char and gas phases along the bed. Nevertheless, for wood gasification, these results also suggest that the main difference between the simulations realized using the Syamlal–O’Brien correlation and those using the Gidaspow and BVK correlations occurs for the predicted length of the particle bed.

Figure 11 shows the axial profiles for the volumetric fraction of the wood phase along the reactor. One may note from viewing Figure 11 that there are almost no differences between the axial profiles for the wood volumetric fraction predicted in the simulations realized with the three drag correlations investigated.

Table 11 presents the simulation results using three different gas–solid heat transfer coefficient correlations. One may note from viewing Table 11 that the RMS values for all simulations are very similar: 0.9% for the Gunn correlation, 0.8% for the Ranz–Marshall correlation, and 0.7% for the Li–Mason correlation. This trend suggests that the gas–solid heat transfer coefficient model has little influence on simulating wood gasification. Concerning the gas temperature leaving the reactor, results in Table 11 show a smaller relative difference when the Ranz–Marshall correlation (Dif_R = 1.7%) is used. The results obtained, when the Gunn and Li–Mason correlations are used, display similar values for the gas temperature leaving the reactor, leading to relative differences of 7.9% and 7.3%, respectively.

Figure 12 presents the axial profiles for the gas temperature along the reactor predicted by the simulations using the three gas–solid heat transfer coefficient correlations, as well as the experimental data from reference [29]. The axial profiles are seen to be in satisfactory agreement with the experimental data in the bed region. Nevertheless, the Ranz–Marshall correlation leads to the best agreement with the experimental data (Figure 12). Furthermore, it is verified that the differences between the gas temperature profiles that originated from the investigated correlations are in agreement with those reported by Xiong and Kong [5] for biomass pyrolysis in bubbling fluidized beds.

Concerning the axial profiles for the volume fractions predicted using different correlations for the gas–solid heat transfer coefficient, only the volume fraction for wood leads to noticeable differences, as can be seen in Figure 13. The axial profiles for the volume fractions of the other solid phases are not presented here since no significant difference is identified. Figure 13 shows that the wood volume fraction predicted using the Gunn correlation, which is smaller than those evaluated in the simulations adopting the Ranz–Marshall and the Li–Mason correlations. This fact is justified because the heat exchange between the gas and wood phases is larger in the simulation using the Gunn correlation than in the simulations adopting the other correlations. Consequently, the temperature of the wood phase calculated in the simulation using the Gunn correlation is higher, as can be seen in Figure 14. Since the wood temperature predicted is higher, wood pyrolysis reaction rates are also higher in this case, leading to a fast consumption of wood. Furthermore, Figure 15 shows that the axial profiles for the char temperature predicted using the Gunn correlation also has higher temperatures than those predicted using the other heat transfer correlations studied here.

It is worth mentioning that, most likely, by changing the formation of the gaseous species and using different choices of chemical mechanisms would presumably alter the influence of the drag and heat transfer models on the simulations’ predictive capabilities. However, this type of investigation is beyond the scope of this work.

4.3. Case 3

Figure 16 and Figure 17 show the biomass volume fraction and the gas temperature profiles along the reactor for the three drag models and three heat transfer correlations that were tested. From Figure 16 and Figure 17, one can note that the drag model choice has more influence on these two variables than the heat transfer correlation.

Figure 16a–c show a segregation of the sugarcane bagasse particles, with the biomass concentrated in the end region of the bed. While the Syamlal O’Brien model predicts a severe segregation, with almost all particles concentrated at the end of the bed, the BVK drag model yields a more uniform biomass distribution along the bed. However, for the BVK drag model, different biomass volume fraction profiles are obtained for the three heat transfer correlations. The greatest difference, as shown in Figure 16c, occurs when the Gunn correlation is used together with the BVK model. In this situation, it cannot be predicted whether the sugarcane bagasse particles will segregate at the end of the bed. The Gunn heat transfer coefficient correlation may have intensified the thermal decomposition of the sugarcane bagasse particles, increasing their consumption rate. This trend is also seen in Case 2 for wood gasification. Conversely, the axial profiles for the biomass volume fraction obtained with the Gidaspow and Syamlal O’Brien drag models seems to be negligibly influenced by the heat transfer correlation. The Gidaspow model leads to a concentration of biomass at the end of the bed, which is higher than the one obtained with the Syamlal–O’Brien model. Consequently, the Gidaspow model predicts less segregation than the Syamlal–O’Brien model.

It is also interesting to note that the different axial profiles in Figure 16a–c reveal a significant impact of the drag correlation on the predicted biomass volume fraction. This finding is also reported by Xiong and Kong [5], who investigated the drag and heat transfer coefficient correlation for a biomass pyrolysis process in a bubbling fluidized bed. However, it contrasts with the findings of Case 2, in which a negligible difference was found (Figure 11). This suggests that the influence of the gas–solid drag correlation on the predicted volume fraction of solid phases is more significant in biomass gasification than in biomass pyrolysis. This trend can be traced to a higher influence of chemical reactions on the granular flow, and also to the corresponding higher temperatures observed in gasification reactors.

Figure 17a–c show the temperature profiles along the reactor height for different drag models and heat transfer correlations. An observation of Figure 17a–c reveals that the BVK drag model used with the Gunn heat transfer correlation provides the highest temperatures along the entire reactor. This behavior is expected since there is a smoother biomass distribution profile along the bed, as shown in Figure 16c. The smoother profile improves the mixture between the phases, and raises the gas temperature. Furthermore, again, due to the more uniform biomass volume fraction profile, the BVK drag model leads to a temperature of the bed that is close to the temperature of the freeboard. On the other hand, the Syamlal–O’Brien drag model leads to a significant drop in the gas temperature from the bed to the freeboard, especially when the Li–Mason correlation is used. This behavior is a consequence of the increased biomass concentration in the end region of the bed. Thus, in this region, the pyrolysis reactions are more intense, causing a sudden drop in the gas temperature.

In

Table 12. Distribution of pyrolysis products for different drag models and heat transfer correlations.

		YNCG [%]	Ytar [%]	Ychar [%]
Experimental [47]		21.00–29.00	62.00–68.00	8.00–10.00
Gidaspow	Gunn	27.36	63.86	8.78
	Ranz–Marshall	27.48	63.86	8.66
	Li–Mason	27.09	64.07	8.85
Syamlal–O’Brien	Gunn	26.51	64.75	8.74
	Ranz–Marshall	26.19	64.91	8.89
	Li–Mason	26.16	64.79	9.05
BVK	Gunn	28.72	62.66	8.62
	Ranz–Marshall	26.82	64.27	8.91
	Li–Mason	27.32	63.85	8.83

and

Table 13. Elementary composition of tar for different drag models and heat transfer correlations.

		C [%]	H [%]	O [%]
Experimental [47]		43.29–47.29	6.30–6.50	46.31–50.31
Gidaspow	Gunn	47.38	6.90	45.73
	Ranz–Marshall	47.33	6.90	45.77
	Li–Mason	47.33	6.92	45.75
Syamlal–O’Brien	Gunn	47.26	6.90	45.84
	Ranz–Marshall	47.21	6.93	45.86
	Li–Mason	47.15	6.94	45.91
BVK	Gunn	47.37	6.90	45.77
	Ranz–Marshall	47.42	6.90	45.63
	Li–Mason	47.48	6.92	45.78

, the distribution of pyrolysis products and the calculated elementary composition of tar at the reactor’s outlet are presented, respectively. Table 12 shows that the drag and heat transfer correlations only have a small impact on the distribution of pyrolysis products. The highest relative difference between the predicted results shown in Table 12 is smaller than 10%, and occurs when there is a sufficient percentage of non-condensable gases present (Y_NCG); this happens when using the Syamlal O’Brien and Li–Mason correlations, and when using the BVK and Gunn correlations. It is possible to observe that the choice of the heat transfer correlation has a greater effect on the distribution of products when the BVK drag model is used. Using the BVK, the highest relative difference found for Y_NCG is approximately 7%, and occurs between the values obtained with the Gunn and the Ranz–Marshall correlations. The variations yielded with the heat transfer correlation for the other drag models are smaller than 2.5%.

Table 13 shows that the elementary composition of tar is also not significantly affected by the choice of drag model and heat transfer correlation. It is interesting to observe that even with great changes in local parameters, such as in the biomass volume fraction and gas temperature, the model is still capable of predicting the distribution of pyrolysis products (Table 12) and the elementary composition of tar (Table 13). In this case, the results suggest that the choice of the interface model or models is only important to evaluate the local parameters of the flow. Unfortunately, there seems to be no experimental information about the distribution of sugarcane bagasse particles and gas temperatures along the reactor, making it difficult to ascertain which closure model yields more accurate axial profiles.

5. Conclusions

The influence of gas–solid drag and heat transfer coefficient models on the prediction capacity of the Euler–Euler approach to simulate reactive and polydisperse bubbling fluidized bed flows is studied. Three different cases are considered, a non-reactive bidisperse bubbling fluidized bed flow (Case 1), and two reactive polydisperse flows in bubbling fluidized beds, one for biomass gasification (Case 2) and the other for biomass pyrolysis (Case 3). The Gidaspow, Syamlal–O’Brien, and BVK gas–solid drag models and the Gunn, Ranz–Marshall, and Li–Mason gas–solid heat transfer correlations are investigated. The computational results for Cases 1, 2, and 3 outline, respectively, the solid mass fraction distribution along the bed, axial profile for gas temperature and the yielded gas composition, and pyrolysis products. All simulation results are compared to appropriate experimental data from the literature.

The results obtained for Case 1 show that the Syamlal–O’Brien and Gidaspow drag models lead to similar results for the non-reactive bidisperse flow, while the BVK drag model yields an axial profile for solid mass fraction which is very different from the one reported in the experimental data. For the reactive biomass gasification process (Case 2), the three drag correlations resulted in similar solid distributions along the reactor. The Syamlal–O’Brien model, however, predicted a smaller bed height than the other models. Furthermore, the Gidaspow model resulted in a more satisfactory prediction of the yielded gas composition and temperature leaving the reactor. Finally, for Case 3, despite the differences in local variables, the predicted compositions of pyrolysis products leaving the reactor using different gas–solid drag and heat transfer coefficient models are very similar. The same is observed for the elementary composition of tar. It is found that the distribution of products is more sensitive to the choice of the drag model when the heat transfer correlation of Gunn is employed, and it is also more sensitive to the choice of the heat transfer correlation when the BVK drag model is employed.

Based on the results presented and discussed, it is possible to conclude that the kinetic modeling of the reactions has a global effect on the simulation of reactive flows in bubbling fluidized bed reactors, which superposes with the influence of the drag and heat transfer coefficient models. Nevertheless, local parameters can be noticeably affected by the choice of the interface closure models.

Finally, based on this study, the Gidaspow drag model and the Ranz–Marshall heat transfer coefficient model can be appropriately adopted in order to simulate reactive flows in bubbling fluidized bed reactors.