A Numerical Study of Turbulent Combustion of a Lignocellulosic Gas Mixture in an Updraft Fixed Bed Reactor

Abstract

Lignocellulosic biomass is an established source of energy with various applications. Yet, its diversity renders the proper combustion of its thermochemical degradation vapors challenging. In this work, the combustion of syngas obtained from biomass thermochemical conversion was numerically investigated to limit pollutant emission. The Computational Fluid Dynamics (CFD) simulation was performed using the open-source OpenFOAM. The reactor was considered in an axisymmetric configuration. The gas mixture resulting from the pyro-gasification devolatilization was composed of seven species: CO, CO₂, H₂O, N₂, O₂, light, and heavy hydrocarbon, represented by methane (CH₄) and benzene (C₆H₆), respectively. The evolutions of mass, momentum, energy, and species’ concentrations were tracked. The flow was modeled using the RANS formulation. For the chemistry, reduced kinetic schemes of three and four steps were tested. Moreover, the Eddy Dissipation Concept (EDC) model was used to account for the turbulence–chemistry interaction. The numerical prediction enabled us to describe the temperature and the species. Results show that all transported variables were closely dependent on the mass flow rate of the inflow gas, the primary and the secondary air injections. Finally, from a process perspective, the importance of the secondary air inlet to limit pollutants emissions can be concluded.

1. Introduction

Crude oil, natural gas, and coal have been used for decades as the primary energy sources. However, excessive use of these fossil fuels has become the main cause of harmful environmental and health effects. Indeed, enormous amounts of greenhouse gas (GHG) are released [1], leading to global warming and brutal climate change. Hence, solar, wind, geothermal, sea waves, and biomass are investigated as viable alternatives and environmentally-friendly renewable energy sources. In particular, biomass has gained increased interest in the last decade due to its availability, abundance, low cost, viability, and sustainability. More precisely, biomass combustion [2,3,4], pyrolysis [5,6,7], and gasification [8,9] have been proven as three processes with high efficiency during energy conversion [10]. In this scope, this work proposes the use of an open-source numerical tool helping in the energy valorization of unconventional woody biomass, such as olive mill solid waste (a substrate much encountered in Northern Africa). Pyrolysis is a widely used thermochemical process. It requires a range of temperatures between 300 °C and 1000 °C and a residence time varying between one to three hours. Different fractions of volatile organic compounds (VOC), bio-oil, tars, and char (fixed carbon and ash) are usually yielded. As for gasification, carbonaceous feedstock should react at high temperatures (higher than 600 °C) in a weak oxidant atmosphere, such as; N₂/O₂ (air), N₂/CO_2, and N₂/steam, depending on the chosen operating conditions. The main product targeted in this case is the syngas composed of CO, H_2, and C_nH_m for which methane is in the majority [11,12]. This gas could be used directly for feeding furnaces, turbines, specific engines [13,14,15], power generation, or for chemical compounds, such as bio-methanol or bio-ethanol [16,17,18]. From an operational perspective, both pyrolysis and gasification could be implemented within the same reactor in a pyro-gasification continuous process. After tars filtration, the rest of the gas could be directly injected into combustion chambers. The gasification process is commonly realized at different isothermal temperature ranges, but always higher than 600 °C and sometimes reaching temperatures as high as 1400 °C. It is to be highlighted that parameters, such as the reactor design, the heating rate, the temperature, the residence time, the gasifier, and the feedstock nature strongly affect the content rates of yielded products [19,20,21,22,23,24,25,26]. Moreover, many reactors could be used, such as the fixed bed, the entrained bed, and the fluidized bed reactors. Fixed bed reactors are reported in the literature in two varieties according to the flow type: the updraft and the downdraft fixed bed configurations [25]. The choice of the adequate gasifier depends on the size requirement (volume of feedstock and energy demand), the feedstock characteristics (moisture and ash content), and the quality of the gasification products (syngas, bio-oils, and bio-char). For example, updraft flow technology is preferred when the feedstock has a moisture content of less than 50%, as in the present study. Moreover, materials with high ash content are not advisable for fixed bed reactors, and woody biomass is not recommended for entrained flow bed reactors [26]. In industrial applications, different feedstocks could be used during the gasification process. However, the most commonly used feedstocks are coal, petroleum-based materials, coke, bitumen, and organic material, such as lignocellulosic biomass and wastes [21,27,28]. Other alternative feedstocks include industrial sludge, fluff, vehicle tires, dried sewage sludge, low-rank coal, and sugarcane bagasse [12]. The present study investigates the combustion of a woody gas mixture similar to that produced by a woody biomass’s pyro-gasification obtained by Farokhi et al. [29] (similar to olive mill solid waste syngas). The combustion of such gas mixtures showed a reduction of NOx emissions during the combustion of CO and H₂ with air [30]. Furthermore, this gas mixture was recently investigated as an alternative fuel for reciprocating engines and gas turbines in the context of heat and/or electricity generation or co-generation processes [31,32]. The Computations Fluid Dynamics (CFD), which is low cost compared to experimental tests, has become a crucial tool for understanding the combustion process of the gas mixture and for controlling the emission levels. In this context, Zhou et al. [33] recently built a model allowing for the reproduction of the combustion of blended solid fuels. The major advantage they concluded was the flexibility of such a model in understanding and predicting the chemical kinetics of reactions. Indeed, after a lab-scale validation, these authors were able to extend the modeling to the scale of an industrial boiler fed with waste wood chips characterized by different ash contents. Inspired by these previous works [29,33], we selected the geometry of an updraft fixed bed reactor where the gas mixture will be burnt above the bed (in the freeboard zone). In a practical case, this configuration avoids the high cost of syngas cleaning from other undesirable components, such as tars.

The primary focus of this study is the effect of the inlet airflow (primary/secondary) and the fuel flow rate on the temperature and species’ concentrations along the axis of the reactor, as that is where experimental measurements were performed. Therefore, a model was built based on Open-FOAM software. The Realizable k-ε model, proven in many studies, was chosen for the turbulence description. Seven species (CO, CO₂, H₂, H₂O, CH₄, and C₆H₆) were considered [33,34]. The Eddy Dissipation Concept model (EDC) was selected with reduced kinetic schemes to consider the turbulence–chemistry interactions. The proposed work would improve the actual reactor’s performance and understand the impact of feedstock variations and the effects of the air supply (primary and secondary) on thermal performances and pollutant emission.

2. Validation Experiments and Feedstock Variability

The experimental data used to validate our calculations correspond to the case of wood pellets used as feedstock in an updraft fixed bed gasifier [29]. In the experimental device, the grate was designed to allow air injection through 32 circular holes (Figure 1a—left, see Section 3.2 for accurate description). After ensuring the beginning of pellet ignition, the decomposition process starts, and the products are ready to react with primary air. The secondary air is injected via 12 holes placed above the grate in an axis-symmetric configuration around the axis. What interests us first of all during calculations is the flammable gas mixture without hydrogen. According to Shen et al. [35], woody biomass is characterized by low ash content (less than 2%). Despite that, it could be higher in some cases, such as for some Nigerian species, Entada Giga 5.09%, and DelonixRegia 4.98% [36]. Woody biomass also presents a high volatile matter content reaching 81.5% (wet basis) for the maritime pine when the moisture content is 8.6% [37]. These properties make woody biomass a high-quality fuel and a good candidate for the pyro-gasification process. It is equally important to underscore that the conversion of woody biomass releases harmful emissions, such as some volatile organic compounds (phenolic aliphatic hydrocarbons), NO_x (NO + NO₂), SO_x (SO + SO₂), and particulate matter emissions (PMs). These emissions are undesirable and must be mitigated as much as possible. It was found, for example, that high amounts of gaseous and particulate emissions were generated by the combustion of the wood in a household boiler [38]. Moreover, hazardous air pollutants (HAPs), depending mainly on the biomass properties, could be deduced after characterization via ultimate analysis allowing the knowledge of the contents of carbon, oxygen, hydrogen, nitrogen, and sulfur. In addition, the proximate analysis provides other crucial properties, such as porosity, density, particle size, fixed carbon content, ash content, moisture content, and the high heating value (HHV) [39]. Additionally, the HAPs are also influenced by the reactor’s design and the operating conditions [40]. Nevertheless, woody biomass presents the lowest nitrogen and sulfur contents among other biomasses and coals, knowing that these elements are directly responsible for NO_x and SO_x emissions [41]. Furthermore, it was reported that the HHV of the gas mixture obtained by gasification depends on the fuel and the gasifier [42,43,44,45]. Furthermore, the low heating value (LHV) is strongly influenced by the operating conditions; the equivalence ratio [42,43], the gasifier temperature [44], the residence time [45], and the moisture content [46]. For information, the LHV value was found to be equal to 13.3 MJ/m³ during the gasification of the wood granulate when using water steam as gasifying agent [44]. According to Couto et al. [47], after analyzing different wood species, the synthesis gas from wood was distinguished by its HHV reaching 14.68 MJ/m³. This diversity of feedstock and resulting properties highlights the need for a free/open-source numerical tool to help engineers to better design/adapt their furnaces to their substrate specificity.

3. Modeling Tools

The combustion of the woody syngas was simulated using the open-source CFD platform Open-FOAM (version 18.12+) [48]. In our case, the study was restricted to only combustion in the gas phase (in the freeboard of the reactor) because we targeted to use this process to feed a combustion chamber directly with this gas mixture [49,50,51].

3.1. Flow, Turbulence, and Combustion Model

RANS equations were used to describe the turbulent flow and combustion of woody biomass syngas. They are composed of the continuity, momentum, energy, and species transport equations [52]. In the current work, the main assumptions are:

• The governing equations are implemented in an unsteady state
• The fuel gas mixture is compressible
• Soret (diffusion) effect, gravitational acceleration effect, and radiant effect are neglected.
• The pressure follows the equation-of-state of an ideal gas (Equation (7))

$$\frac{\partial \overline{\rho }}{\partial t}+{\displaystyle \mathsf{\nabla }}.\left(\overline{\rho }\tilde{U}\right)=0$$

(1)

$$\frac{\partial \left(\overline{\rho }\tilde{U}\right)}{\partial t}+{\displaystyle \mathsf{\nabla }}.\left(\overline{\rho }\tilde{U}\otimes \tilde{U}\right)-{\displaystyle \mathsf{\nabla }}.\tau =-{\displaystyle \mathsf{\nabla }}\overline{p}$$

(2)

$$\frac{\partial \left(\overline{\rho }\tilde{h}\right)}{\partial t}+{\displaystyle \mathsf{\nabla }}.\left(\overline{\rho }\tilde{U}\tilde{h}\right)+\frac{\partial \left(\overline{\rho }\tilde{k}\right)}{\partial t}+{\displaystyle \mathsf{\nabla }}.\left(\overline{\rho }\tilde{U}\tilde{k}\right)-\frac{\partial \overline{p}}{\partial t}={\displaystyle \mathsf{\nabla }}.\left({\alpha }_{eff}{\displaystyle \mathsf{\nabla }}\tilde{h}\right)+{\overline{Q}}_{\mathit{heat}}$$

(3)

$$\frac{\partial \overline{\rho }\widetilde{Y_i}}{\partial t}+{\displaystyle \mathsf{\nabla }}.\left(\overline{\rho }\tilde{U}\widetilde{Y_i}\right)={\displaystyle \mathsf{\nabla }}.\left({\mu }_{eff}{\displaystyle \mathsf{\nabla }}\widetilde{Y_i}\right)+\overline{R_i}$$

(4)

$$\overline{p}=\overline{\rho }RT {{\displaystyle \sum }}_i\frac{\widetilde{Y_i}}{\overline{w_i}}$$

(5)

$${\overline{Q}}_{\mathit{heat}}=-{\displaystyle \mathsf{\nabla }}.\left(\lambda {\displaystyle \mathsf{\nabla }}T\right)+\overline{\rho }{{\displaystyle \sum }}_{i=1}^N\widetilde{h_i}\widetilde{Y_i}\widetilde{U_i}$$

(6)

$$\tau =\overline{\rho }{\mu }_{eff} {\displaystyle \mathsf{\nabla }}\tilde{U}+\overline{\rho }{\mu }_{eff}\left[{\left({\displaystyle \mathsf{\nabla }}\tilde{U}\right)}^{T_r}-\frac{2}{3}\left({\displaystyle \mathsf{\nabla }}.\tilde{U}\right)I\right]$$

(7)

where $\overline{\rho }$ is the mean density, $\tilde{U}$ is the Favre-averaged velocity, t is the time, $\tau$ is the averaged stress tensor, $\overline{p}$ is the mean pressure, $\tilde{h}$ is the Favre-averaged internal energy, $\tilde{k}$ is the Favre-averaged specific kinetic energy, ${\alpha }_{eff}$ is the effective thermal diffusivity, ${\overline{Q}}_{heat}$ is the mean heat generation rate, μ_eff is the effective dynamic eddy viscosity, R is the gas constant, T is the temperature, $\overline{R_i}$ is the mean reaction rate, $\widetilde{Y_i}$ and $\overline{w_i}$ are the Favre-averaged mass fraction and mean molar mass for the ith species, λ is the thermal conductivity of the mixture, T_r refers to the transpose of ${\displaystyle \mathsf{\nabla }}\tilde{U}$ and I is the identity tensor. Turbulence was modeled using the Realizable k-ε turbulence model [53,54,55]. The transport equations for the turbulent kinetic energy k and its dissipation rate ε in this model are written as [56]:

$$\frac{\partial \overline{\rho }\tilde{k}}{\partial t}+{\displaystyle \mathsf{\nabla }}.\left(\tilde{U}\overline{\rho }\tilde{k}\right)-{{\displaystyle \mathsf{\nabla }}}^2\left(\overline{\rho } D_{k_{eff}}\tilde{k}\right)=\overline{\rho }G-\frac{2}{3}\overline{\rho }\left({\displaystyle \mathsf{\nabla }}. \tilde{U}\right)\tilde{k}-\overline{\rho }\epsilon +S_k$$

(8)

$$\frac{\partial \overline{\rho }\tilde{\epsilon }}{\partial t}+{\displaystyle \mathsf{\nabla }}.\left(\tilde{U}\overline{\rho }\tilde{\epsilon }\right)-{{\displaystyle \mathsf{\nabla }}}^2\left(\overline{\rho }{D}_{{\epsilon }_{eff}}\tilde{\epsilon }\right)=C_1\overline{\rho }\left|S\right|\tilde{\epsilon }+C_2\overline{\rho }\frac{{\tilde{\epsilon }}^2}{k+\sqrt{\nu \tilde{\epsilon }}}+S_{\epsilon }$$

(9)

where $\tilde{k}$ represents the Favre-averaged of turbulent kinetic energy, $\tilde{\epsilon }$ represents the Favre-averaged turbulent dissipation rate, G is the turbulent kinetic energy production rate due to the anisotropic part of the Reynolds-stress tensor, D_εeff is the effective diffusivity for ε, D_keff is the effective diffusivity for k. C₁ is the model coefficient (C₁ =1.44), C₂ is the model coefficient equal to 1.92, S_k is the internal source term for k, S_ε is the internal source term for ε, and S is the mean strain rate. The expression of the boundary conditions for the turbulent kinetic energy k and its dissipation rate ε follow the expressions below [49]:

$$k=\frac{3}{2}{\left(U\cdot I\right)}^2$$

(10)

$$\epsilon =\frac{C_{\mu }^{0.75}\cdot k^{\frac{3}{2}}}{L_t}$$

(11)

where, U and I are the flow velocity and the turbulence intensity, respectively. C_µ and L_t are the turbulent intensity and its integral length scale.

The EDC approach [57,58] was proven a potent tool for investigating combustion in moderate or low oxygen dilution conditions. The idea consists in dividing the computational domain into reacting fine structures cells and non-reacting regions called surroundings. The reacting regions are modeled as perfectly stirred reactors (PSR) needing fine structure treatment since reactants are mixed on a molecular scale, and mass is exchanged with the surroundings. The fine structures are the seats of the smallest scales of turbulence. If turbulence and turbulent mixing decrease, the reaction progress can be limited by either mixing or chemistry. In such a regime, finite rate chemistry is necessary to accurately describe the turbulence-chemistry interaction [59]. However, those reactions could extend outside the fine structures due to the species’ mass transfer. Therefore, the fine structure mass share is expressed as [60]:

$${\gamma }^{*}=C_{\gamma }^{\frac{1}{2}}{\left(\frac{\nu \epsilon }{k^2}\right)}^{\frac{1}{2}}={\gamma }_L^2$$

(12)

where, $C_{\gamma }={\left(\frac{3C_{D2}}{4C_{D1}^2}\right)}^{1/4}\approx 2.13$. C_D1 equal to 0.135 and C_D2 equal to 0.5 are two constants of the EDC model [61]. $\nu$ is the turbulent viscosity. $\epsilon$ is the turbulence dissipation rate, and k is the turbulence kinetic energy. As proposed by Magnussen [60], the mean reaction rate between a certain fraction χ of fine structure regions and the rest of the fluid ($\overline{R_i}$) is given by:

$$\overline{R_i}=\frac{\overline{\rho }}{{\tau }^{*}}\frac{{\gamma }_L^2\chi }{1-{\gamma }_L^2\chi }\left(\widetilde{Y_i}-Y_i^{*}\right)$$

(13)

Whereas, the fine structure and the surroundings are linked by the fluid averaged properties and the fine structure mass share ${\gamma }^{*}$ as:

$$\overline{\psi }={\psi }^0\left(1-{\gamma }^{*}\right)+{\psi }^{*}{\gamma }^{*}$$

(14)

In the above Equations (14)–(16); *, ◦, − and $~$ denote, respectively, the fine structure, the surrounding, and the fluid mean properties (Reynolds and Favre averaged), Y_i* is the fine structure mass fractions of species i. χ represents the reacting fraction of the fine structures. τ* is the fine structures residence time or EDC mixing time. γ_L is the fine structure length fraction. γ* is the fine structure share.

According to the original EDC’s version [57], the fine structures residence time τ* and a dimensionless fine structure length fraction could be defined based on the Reynolds number as follows:

$$R{e}_t=\frac{k^2}{\epsilon \nu }$$

(15)

$${\gamma }_L=C_{\gamma }R{e}_t^{-1/4}$$

(16)

$${\tau }^{*}=C_{\tau }R{e}_t^{-1/2}\frac{k}{\epsilon }$$

(17)

With, $C_{\tau }={\left(\frac{C_{D2}}{3}\right)}^{1/2}.$

To manage the computational cost of integrating the chemistry, only reduced kinetic schemes could be considered. All these equations were implemented in the reacting-Foam solver (part of OpenFOAM).

3.2. Simulation Setup

The computational domain is an updraft fixed bed reactor in the freeboard zone, operating at a nominal capacity of 8–11 kW. The combustion chamber is a cylinder with a diameter of 200 mm and a height of 820 mm, as shown in Figure 1a. Three inlets characterize this reactor: the first is reserved for the woody biomass syngas feed and is composed of 32 holes of 3 mm diameter. The second and the third serve to supply the reactor with primary and secondary air with 12 holes of 4 mm diameter. These two air inlets are placed above the bed at 60 mm and 75 mm, respectively. Figure 1b presents a sketch of the down part of the reactor bed with a 2D wedge mesh. More details of the reactor design were reported in the literature [62,63]. Such setups are notoriously greedy computational power, which would hinder engineers’ reuse of our work. Therefore, the reactor geometry was simplified by adopting an axisymmetric 2D configuration [52]. Therefore, the holes of the primary and secondary air inlets are not resolved. This emulates a concentrically slid inlet, which, as it has a different inlet area, leads to smaller momentum of inflows. To remedy this shortcoming, a non-uniform grid mesh was adopted with a fine mesh near the holes of primary and secondary air. The grid mesh geometry’s general view is shown in Figure 1b, where a non-uniform mesh is used. Thus, the computational domain contains 32,772 (see Section 3.3 for convergence analysis).

Once the model is set and the geometry meshed, the simulation setup (boundary conditions, chemical reaction parameters…) is to be considered. In this regard, it is to be highlighted that it is difficult to reproduce precisely the experimental conditions because of uncertainties in the flow rate values and values of k and ε. Three cases named Case 1, Case 2, and Case 3 were investigated to identify which gives the best agreement with experimental results obtained by Farokhi et al. [29]. Data characterizing these cases are summarized in

Table 1. Mass flow rate, turbulence kinetic energy k (m²·s⁻²), and its dissipation rate ε (m²·s⁻³). The initial pressure condition for gas and air is P = 101,325 Pa. P. air—Primary Air, S. air—Secondary Air.

	Case 1			Case 2			Case 3
$\dot{m} \left(\mathrm{Gas}\right)\left({10}^{-6} \mathrm{kg}·{\mathrm{s}}^{-1}\right)$	6.0			5.0			6.0
$\dot{m} \left(\mathrm{P}.\mathrm{air}\right)\left({10}^{-6} \mathrm{kg}·{\mathrm{s}}^{-1}\right)$	0.20			0.20			0.16
$\dot{m} \left(\mathrm{S}.\mathrm{air}\right)\left({10}^{-6} \mathrm{kg}·{\mathrm{s}}^{-1}\right)$	1.50			1.0			2.50
Kinetic energy of turbulence k$({\mathrm{m}}^2·{\mathrm{s}}^{-2})$	Gas 3.18 × 10−2	P. air 2.98 × 10−5	S. air 1.06 × 10−5	Gas 2.21 × 10−2	P. air 2.98 × 10−5	S. air 0.47 × 10−5	Gas 3.18 × 10−2	P. air 1.90 × 10−5	S. air 2.96 × 10−5
Dissipation rate of turbulence $\epsilon \left({\mathrm{m}}^2·{\mathrm{s}}^{-3}\right)$	Gas 1.33 × 10−2	P. air 6.68 × 10−5	S. air 4.51 × 10−5	Gas 0.77 × 10−2	P. air 6.68 × 10−5	S. air 1.33 × 10−5	Gas 1.33 × 10−2	P. air 3.40 × 10−5	S. air 20.92 × 10−5
Flow velocity U (m·s−1)	Gas 2.90	P. air 0.61	S. air 1.45	Gas 2.41	P. air 0.61	S. air 0.97	Gas 2.9	P. air 0.49	S. air 2.42

. In order to have an idea of the fine structures residence time for each case simulation, Equations (17) and (19) were used to calculate this crucial parameter a priori (primary air and secondary air). Indeed, the reactor residence time (L/U_max ≈ 0.3 s) should be much larger than the fine structure residence time to ensure the validity of the EDC turbulence modeling approach.

Table 2. Calculation of turbulent Reynolds number and the fine structures residence time for each case.

Gas Mixture				Primary Air				Secondary Air
Case 1				Case 1				Case 1
K 3.18	$\epsilon$ 1.33	Ret 3.80 × 106	${\tau }^{*}$ 5.0 × 10−4	K 2.98	$\epsilon$ 6.68	Ret 6.6 × 105	${\tau }^{*}$ 5.4 × 10−4	K 1.06	$\epsilon$ 4.51	Ret 1.2 × 105	${\tau }^{*}$ 2.7 × 10−4
Case 2				Case 2				Case 2
K 2.21	$\epsilon$ 0.77	Ret 3.17 × 106	${\tau }^{*}$ 6.5 × 10−4	K 2.98	$\epsilon$ 6.68	Ret 6.6 × 105	${\tau }^{*}$ 5.4 × 10−4	K 0.47	$\epsilon$ 1.33	Ret 0.83 × 105	${\tau }^{*}$ 5.0 × 10−4
Case 3				Case 3				Case 3
K 3.18	$\epsilon$ 1.33	Ret 3.8 × 105	${\tau }^{*}$ 5.0 × 10−4	K 1.90	$\epsilon$ 3.40	Ret 5.3 × 105	${\tau }^{*}$ 3.1 × 10−4	K 2.96	$\epsilon$ 20.92	Ret 2.1 × 105	${\tau }^{*}$ 1.2 × 10−4

exhibits the values of residence time corresponding to the three cases and based on the gas mixture properties presented in Table 1.

It could be concluded from Table 2 that the residence time of secondary air was the lowest in Case 1 and Case 2 (2.7 × 10⁻⁴ s and 1.2 × 10⁻⁴ s, respectively). For the rest of the cases, residence times were in the vicinity of 5 × 10⁻⁴ s. All of them are much lower than the reactor residence time, allowing us to conclude the validity of the EDC approach for the three flow rate configurations. Moreover, on top of uncertainty regarding the experimental setup, the reaction schemes used in these calculations are also a source of concern. Therefore, different chemical mechanisms were considered to evaluate the temperature and gaseous emissions to determine the most adequate.

Table 3. Chemical kinetic mechanisms and corresponding Arrhenius coefficients.

Mechanism Number	Chemical Kinetics Mechanism	A (s−1)	β	Ta (K)	Reaction Order	References
T-S.1	CH4 + 1.5O2 = CO + 2H20	1.398 × 10+10	−0.062	14,038	[CH4]0.5 [O2]1.066	[64]
	CO + 0.5O2 = CO2	7.66 × 10+11	0.215	9213	[CO] [O2]	[64]
	C6H6 + 4.5O2 = 6CO + 3H2O	2.4 × 10+11	0	15,098	[C6H6]−0.1 [O2]1.85	[65]
T-S.2	CH4 + 1.5O2 = CO + 2H2O	1.400 × 10+10	−0.062	14,040	[CH4]0.5 [O2]1.066	[31]
	CO + 0.5O2 = CO2	7.380 × 10+11	0.215	9209	[CO] [O2]0.5	[31]
	C6H6 + 4.5O2 = 6CO2 + 3H2O	1.7 × 10+11	0	15,098	[C6H6]−0.1 [O2]1.85	[65,66]
F-S.1	CH4 + 1.5O2 = CO + 2H2O	1.4 × 10+10	−0.062	14,038	[CH4]0.5 [O2]1.066	[64]
	CO + 0.5O2 = CO2	2.24 × 10+12	0	20,471	[O2]0.25	[67]
	CO2 = CO + 0.5O2	5 × 10+8	0	20,471	[CO2]	[67]
	C6H6 + 4.5O2 = 6CO2 + 3H2O	2.4 × 10+11	0	15,098	[C6H6]−0.1 [O2]1.85	[65,66]
F-S.2	CH4 + 0.5O2 = CO + 2H2	4.4 × 10+11	0	15,034	[CH4]0.5 [O2]1.25	[66]
	CO + 0.5O2 = CO2	7.66 × 10+11	0.215	9213	[CO] [O2]0.25	[31,68,69]
	H2 + 0.5O2 = H2O	6.8 × 10+15	−1	20,086	[H2]0.25 [O2]1.5	[65,66]
	C6H6 + 4.5O2 = 6CO + 3H2O	2.4 × 10+11	0	15,098	[C6H6]−0.1 [O2]1.85

summarizes the different used chemical kinetic mechanisms. These mechanisms are characterized by the frequency factor (A), the activation temperature (T_a), the temperature exponent (β), and the reaction order. The chemical kinetics were modeled using the Arrhenius law expressed as:

$$\dot{\omega }=A{T}^{\beta }\exp \left(-T_a/T\right){\left[fuel\right]}^b{\left[oxidizer\right]}^c$$

(18)

$\dot{\omega }$ is the reaction rate, b and c are the exponents of the fuel concentration and the oxidizer concentration, respectively.

Neglecting radiative transfer of the fine structure, the chemical species reaction rate ${\dot{\omega }}_i^{*}$ is inherent to the fine structure mass fractions of species $Y_i$ via the following expression:

$$\frac{d{Y}_i}{dt}={\dot{\omega }}_i^{*}+\frac{1}{{\tau }^{*}}\left(Y_i^0-Y_i^{*}\right)$$

(19)

Hence, after resolving Equation (19), the calculation of the mean reaction rate $\overline{R}$ of Equation (15) is then achieved.

All considered mechanisms are reduced scheme types composed of three or four-step global reactions with six or seven species. Three-step mechanism (T-S) was considered for C₆H₆, CH₄ and CO reactivity’s. Four-step (F-S) mechanism was devoted for CO₂ decomposition and was coupled to the three-step mechanisms. All used mechanisms are summarized in Table 3.

Finally, the gas composition is to be addressed.

Table 4. Details of the inlet boundary conditions and gas flow characteristics [29].

Species	C6H6	CH4	CO2	O2	H2O	CO	N2	Primary Air		Secondary Air
Mass fraction	0.130	0.029	0.264	0.019	0.109	0.015	0.432	0.233	0.767	0.233	0.767
Temperature (K)	1373							300		300
Intensity of turbulence (%)	5							23		58
Turbulence length scale L(m)	7 × 10−2							4 × 10−2		4 × 10−2

presents the initial conditions of the woody biomass syngas flow, the two air supply types, the mass fractions of species, the temperature, and the turbulence intensity. The gas mixture properties presented in this table were the results of experimental measurements corresponding to the steam wood pellets pyro-gasification [29].

3.3. Mesh Convergence

Once the numerical setup had been defined and before starting the numerical study, the mesh size was tested to optimize the computation time, while establishing a convergent solution [70]. For this purpose, different (hexahedral) meshes of different sizes were tested until reaching stable results under further mesh refinement. Thus, the meshes named: M1 = 6994 cells, M2 = 17,570 cells, M3 = 32,772 cells, M4 = 45,302 cells, M5 = 58,193 cells, and M6 = 67,010 cells, were tested. Figure 2 shows the flame temperature profile along the axis of the reactor for the different refinements (Case 1 with chemical kinetic mechanism F-S.1, allegedly the most demanding setup). The results of these calculations show that the refinement of the mesh influences the simulation only for the first three types of mesh M1, M2, and M3. Therefore, to save computation time and memory, it is not necessary to further refine the mesh beyond M3. Therefore, the M3 mesh characterized by 32,772 cells must be selected to perform the following calculations.

4. Simulation Results and Discussion

Before diving deep into the results, a first qualitative analysis is proposed. Figure 3a,b show the temperature 2D surface plot in the reactor for Case 1 and Case 3 (the two most different boundary conditions). First of all, they show no sign of numerical artifact. Second, they allow for to identification of the effect of the boundary conditions on the temperature distribution inside the reactor. For example, the oxidation zone (OZ) length is more prominent (approximately +30%) when using Case 1 compared to Case 3. This can be explained by the fact that the mass flow rate, the kinetic turbulence, and the flow velocity U of primary air increase the oxidation zone length along the axis of the reactor. On the contrary, the radial aspect of the combustion zone is more developed in Case 3 than in Case 1 (the radius of OZ-1 is wider than the radius of OZ-3 for T > 800). Indeed, in this case, the larger flow rate in the secondary injection confines the reaction zone.

4.1. Axial Temperature Profile

Figure 4 shows the flame’s temperature variation along the axis of the reactor for both experiments and numerical simulations. The trend is consistent between the simulations and agrees with the experiments. The temperature peak ranges between 1958 K and 2142 K, which is small and was expected given the tested conditions and chemical models [71]. In more detail, decreasing the gas injection is responsible for most of the change in the peak temperature value for a given kinetic mechanism. Manipulation of the secondary air injection has only a marginal effect. Moreover, the choice of the kinetic mechanism also influences this indicator (e.g., a 50 K decrease when switching from T-S.1 to F-S.1, in Case 1). Turbulence level also explains part of the differences observed between the simulations. Indeed, the lower efficiency of the micro-scale mixing between the gas and oxygen can explain the enlargement of the temperature profile at the expense of its sharpness in Case 2 [72,73]. Similar results were reported for other feedstocks combustion in different reactors [74,75,76]. All in all, the best agreement was obtained with the conditions of Case 3 with the F.S-1 kinetic mechanism. Nevertheless, all the tested cases can be considered to represent adequately the experimentally reported temperature.

4.2. Chemical Species Concentrations

Figure 5 illustrates the evolution of CO₂, O₂, and CO concentrations along the axis of the reactor. Even though simulation curves surround the experimental values, these are different from the experimental curves. CO₂ and O₂ predictions are more influenced by the operating conditions (mainly gas rate flows and turbulence characteristics described in each Case), while the reaction mechanism controls CO concentration. As species’ concentrations exhibit significant differences between numerical predictions (while temperature did not), they can be used to assess the most suited model. In our case, as pollutants (e.g., CO) are our focus, it can be concluded that the best results are achieved in Case 1 (best agreement on average). The question of the chemical mechanism remains open. While the T.S-2 kinetic mechanism provides the best agreement on CO final value (length > 0.4 m), its value is under predicted over the whole reactor and drops earlier than experimental results. T.S.-1 slightly overestimates the final CO concentration but provides reasonable agreement over most of the reactor length (length < 0.4 m). As any experiment can show some absolute value measurement error, predicting a good dynamic is to be favored. Therefore, it can be concluded that the T.S.-1 mechanism is the most suited one in this configuration. Diving deeper into the results, the comparison between Case 1 and Case 3 shows that a reduction of the mass flow rate of the primary air in Case 3 leads to less efficient combustion (lower CO₂ concentration and higher O₂ concentration at the same time). This confirms the assumption that reducing the mass flow rate of the primary air (from Case 1 to Case 3) results in a reduction of CO emissions. Moreover, when reducing the mass flow rate of the inflow gas in Case 2, an increase of the mole fraction of CO₂ and the concentration of CO against a large decrease of the O₂ was observed. The increase of the CO₂ indicates that the combustion in Case 2 was more complete and more stable than in Case 1. These results exemplify how CFD can be used to choose operating conditions to lower pollutant emissions.

4.3. Effect of the Mass Flow Rate on the Velocity Field

Figure 6 shows the velocity vectors and the temperature contours for Case 1 when adopting the T-S.1 chemical mechanism. Two recirculation zones, denoted Z₁ and Z₂, were observed. The first zone (Z₁) appears at 0.25 m, where the temperature reaches approximately 1580 K. The interaction between the two supplied air (primary and secondary) in zone Z₁ with the main gas is intensive. An advantageous consequence is the transport of heat toward the burner walls, where it can be recovered. The second zone, Z₂, is situated at the height of 0.35 m, where the temperature reaches 1328 K. In this zone (Z₂), the interaction between the air and the main gas flow appears less intensive, and this gas behavior was reported by Farokhi et al. [29]. More precisely, the latter authors noted that air interactions produce a resistance to the reacting gas, which causes a split of the coming flow into the two up-mentioned zones (Z₁ and Z₂).

4.4. Heat Generation Spatial Distribution

Figure 7 illustrates the mean heat generation rate (${\overline{Q}}_{heat}$) map for Case 1 when adopting T-S.1 mechanism. We note that ${\overline{Q}}_{heat}$ contours are placed at the gas inflow and after the primary air injection. This result correlates well with the temperature profile (Figure 3). At the gas inlet, ${\overline{Q}}_{heat}$ is less than 5000 kJ/kg, whereas it reaches 20,000 kJ/kg after the injection of the primary air. Nevertheless, as one can see, primary air injection is not sufficient to oxidize all the syngas. Indeed, while its flow rate is already high and the combustion reaction intense, it fails to reach the center of the reactor. Therefore, it highlights that proper secondary air injection tuning is mandatory to ensure adequate combustion and oxidation of the supplied gas.

5. Conclusions

This paper presents a full numerical description of the syngas combustion from the water steam wood pellets pyro-gasification in an updraft fixed bed reactor based on the open-source OpenFOAM. The temperature, the gaseous emission evolution of CO, CO₂, and O₂, the velocity field, and the release heat rate were studied. Each simulation run was realized based on a combination of cases (reported in Table 1 and describing the gas rate flows and the turbulence regimes) and a chemical kinetic mechanisms reported in Table 3. It was revealed that:

• Case 3-FS1 was the best combination for predicting the temperature profile despite the peak temperature being under estimated and for all combinations. This may be due to not taking into account of radiative transfer during our calculations;
• Case 1-TS1 was the best for predicting CO₂ emissions;
• Case 1-FS2, Case 1-TS1 and Case 1-FS1 were better for predicting O₂ concentrations;
• Case 1-TS1 predicted better CO concentration for a length < 0.4 m, but Case 2-TS2 was the best one for length > 0.4 m.

Hence, this realized workflow was shown to be capable of reproducing the behavior of an actual burner. Indeed, it was possible to assess the effect of operating conditions (gas inflow, primary and secondary air inflow) and to identify the relevant kinetic model for the tested biomass. Taking a step, this work illustrates how an open-source CFD tool could be used to improve reactor design and operation to increase efficiency and mitigate pollutant generation.