# Simulation of Wood Combustion in PATO Using a Detailed Pyrolysis Model Coupled to fireFoam

Hermes Scandelli, Azita Ahmadi-Senichault, Franck Richard, Jean Lachaud

<figure><img src="image_1.png" /></figure>

## To cite this version:

Hermes Scandelli, Azita Ahmadi-Senichault, Franck Richard, Jean Lachaud. Simulation of Wood Combustion in PATO Using a Detailed Pyrolysis Model Coupled to fireFoam. Applied Sciences, 2021, 11 (22), pp.10570. 10.3390/app112210570. hal-03466603

HAL Id: hal-03466603

https://hal.science/hal-03466603

Submitted on 6 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Simulation of Wood Combustion in PATO Using a Detailed Pyrolysis Model Coupled to fireFoam

Hermes Scandelli<fnref n="1" />*, Azita Ahmadi-Senichault<fnref n="1" />, Franck Richard<fnref n="2" /> and Jean Lachaud<fnref n="1" />

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Scandelli, H.;
Ahmadi-Senichault, A.; Richard, F.;
Lachaud, J. Simulation of Wood
Combustion in PATO Using a
Detailed Pyrolysis Model Coupled to
fireFoam. *Appl. Sci.* **2021**, *11*, 10570.
https://doi.org/10.3390/app112210570

Academic Editor: Maria Grazia De
Giorgi

Received: 14 September 2021
Accepted: 5 November 2021
Published: 10 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

**Abstract:** The numerical simulation of fire propagation requires capturing the coupling between wood pyrolysis, which leads to the production of various gaseous species, and the combustion of these species in the flame, which produces the energy that sustains the pyrolysis process. Experimental and numerical works of the fire community are targeted towards improving the description of the pyrolysis process to better predict the rate of production and the chemical nature of the pyrolysis gases. We know that wood pyrolysis leads to the production of a large variety of chemical species: water, methane, propane, carbon monoxide and dioxide, phenol, cresol, hydrogen, etc. With the idea of being able to capitalize on such developments to study more accurately the physics of fire propagation, we have developed a numerical framework that couples a detailed three-dimensional pyrolysis model and fireFoam. In this article, we illustrate the capability of the simulation tool by treating the combustion of a wood log. Wood is considered to be composed of three phases (cellulose, hemicellulose and lignin), each undergoing parallel degradation processes leading to the production of methane and hydrogen. We chose to simplify the gas mixture for this first proof of concept of the coupling of a multi-species pyrolysis process and a flame. In the flame, we consider two separate finite-rate combustion reactions for methane and hydrogen. The flame evolves during the simulation according to the concentration of the two gaseous species produced from the material. It appears that introducing different pyrolysis species impacts the temperature and behavior of the flame.

**Keywords:** numerical simulations; wood; pyrolysis; combustion; puffing effect

## 1. Introduction

Charring materials (timber and composite materials) are widely used in many civil [1] and aeronautical [2] applications. On the one hand this leads to a sustainable development with lower energy requirements and less pollution, but on the other hand, this may lead to fire hazards. These materials have a flammable nature, and this makes necessary the prediction of a fire potential with the consequent loss of stiffness, strength and resistance of a structure [3]. In the recent few decades, different studies approached the topic and aimed to characterize the performance of charred materials under combustion. Fire retardancy [4], flammability properties [5] and gaseous emissions [6] are just a few examples of the main properties that have been studied. The main strength of this fire community lies in its ability to describe how a fire behaves for a given material, including the description of ignition, fire growing, propagation and extinction processes [7]. Different numerical simulation programs have been developed for such purposes. One of the most successful is fireFoam [8], an opensource code to simulate and visualize the fire behavior in a defined geometry. Different application cases can be found in the literature [9–11].

In a numerical simulation, modeling the flame is not enough. A description of the charred material degradation process in the absence of oxygen (the pyrolysis), as well as the coupling between material and environment, are fundamental aspects that need

<footnote n="1">¹ Arts et Métiers Institute of Technology, University of Bordeaux, CNRS, Bordeaux INP, INRAE, I2M Bordeaux, 33400 Talence, France; azita.ahmadi-senichault@u-bordeaux.fr (A.A.-S.); jean.lachaud@u-bordeaux.fr (J.L.)</footnote>

<footnote n="2">² Institut Pprime (UPR 3346 CNRS), Université de Poitiers, Futuroscope Chasseneuil CEDEX, 86000 Poitiers, France; franck.richard@univ-poitiers.fr</footnote>

<footnote>* Correspondence: hermes.scandelli@u-bordeaux.fr</footnote>
to be captured. Pyrolysis leads to the production of a large variety of chemical species: water, methane, propane, carbon monoxide and dioxide, phenol, cresol, hydrogen, etc. The production of these species and their concentration is strongly affected by the imposed conditions [12,13]. This would not affect only the material, but even the combustion process. The gaseous species, once produced, will percolate in the material and eventually reach the external environment. Many of these species are fuels, therefore the correct prediction of their formation could change the outcome of the simulation. Different pyrolysis models have been developed during the years, each characterized by a different level of accuracy in the description of the chemical reactions. The simplest one is the single reaction mechanism [14] in which only a single decomposition reaction is considered. A next step is done with the multi-component single reaction model [15] where the description becomes sensible to the feed-stock. With an increase in complexity, there is the competitive model [16] that introduces competitive reactions to predict different product distributions depending on the conversion conditions, such as heating rate and pressure. Finally, there are the competitive multi-component models [17,18] which represent efficiently the dependence on the feed-stock as well as the temperature effect on the yields.

Pyrolysis is not the only process to model in charring materials. Mass, momentum and energy of the material need to be conserved, and the transport of gaseous species also need to be taken also into account. A generic combustion problem with a charring material is schematized in Figure 1.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** A generic combustion problem. Two different regions are present: the environment, where combustion happens, and the material, characterized by pyrolysis.</figcaption></figure>

The problem may be split in two regions: the environment, where combustion happens, and the material, characterized by a thermal degradation process, the pyrolysis. Two different models are required for the numerical resolution. We chose to couple fireFoam with a detailed pyrolysis model. The latter enters in the type-3 models classification proposed by Lachaud et al., 2011 [19], and it has already been applied to different applications, like ablative heat-shield design and pyrolysis of lignocellulosic biomass [20–22].

The main objective of this paper is to propose for the first time a numerical framework that couples a detailed three-dimensional pyrolysis model and fireFoam. The numerical framework with the description of the two numerical models is given in the next sections: Section 2 presents the detailed pyrolysis model and Section 3 presents the combustion process. The coupling conditions between the two models are given in Section 4. The capabilities of the solver are showed in Section 5, where two different applications are considered. The first one aims to show how the introduction of different pyrolysis species impacts the temperature and the behavior of the flame. The second application is a 2D simulation of a wood log combustion. Here, we chose to simplify the gas mixture for this first proof of concept of coupling of a multi-species pyrolysis process and a flame. Results show a perfect coupling between the two models and highlight the impact of the species, whose concentration change in time, on the fire behavior. Finally, conclusions are drawn in Section 6.
## 2. Numerical Model: Material Region

The model to describe the material region is presented in this section. It consists of a generic pyrolysis model that allows the description of the interaction between several solid phases and a gas phase [22]. A summary of the governing equations of the solid phases, gas and species are presented in the following subsections.

### 2.1. Main Assumptions

The model considers the interaction between a multi-phase reactive material ($N_p$ solid phases) with a multi-species reactive gas mixture ($N_g$ gaseous elements/species). No liquid phase is modeled. Any liquid present in the material (water) is modeled as a solid phase. The description is done at the macroscopic scale, where the governing equations can be derived from upscaling theories. Their derivation relies on the existence of a representative elementary volume (REV) of the domain and on the assumption of scales separation, as illustrated in Figure 2.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Illustration of the problem: interaction between a multi-phase reactive material with a multi-species reactive gas mixture. The model characterizes this problem at the macroscopic scale under the assumption of scale separation ($L \gg l$).</figcaption></figure>

### 2.2. Pyrolysis

The material is assumed to be composed of $N_p$ solid phases. For example, in the wood-cell walls the main components are cellulose, hemicellulose, and lignin, which are modeled as three phases. Each solid phase $s_i$ may decompose following several pyrolysis kinetics in the absence of oxygen. For this reason, each phase $i$ is split into $P_i$ sub-phases to model the different degradation mechanisms. A generic decomposition of a sub-phase $j$ (from solid phase $i$) leads to the production of species $A_k$ according to the stoichiometric coefficients $\bar{\xi}_{i,j,k}$ as follows

$$s_{i,j} \rightarrow \sum_{k=1}^{N_g} \bar{\xi}_{i,j,k} A_k, \quad \forall i \in [1, N_p], \forall j \in [1, P_i] \quad (1)$$

where $N_g$ is the total number of the gaseous species accounted for in the gas mixture. Pyrolysis reactions are modeled by Arrhenius laws. In this way, the advancement of the pyrolysis reaction $\chi_{i,j}$ of sub-phase $j$ within phase $i$ can be defined as

$$\partial_t \chi_{i,j} = (1 - \chi_{i,j})^{m_{i,j}} T_g^{n_{i,j}} A_{i,j} \exp\left(-\frac{\mathcal{E}_{i,j}}{RT_g}\right), \quad \forall i \in [1, N_p], \forall j \in [1, P_i] \quad (2)$$

where $R$ is the perfect gas constant, $T_g$ the temperature of the gas, $A$ the Arrhenius law pre-exponential factor, $\mathcal{E}$ the Arrhenius law activation energy, and $m$ and $n$ the Arrhenius law parameters. The total advancement in the pyrolysis process, $\tau$, is evaluated as follows
$$ \tau = \sum_{i \in [1, N_p]} \sum_{j \in [1, P_i]} \frac{\epsilon_{i,0} \rho_{i,0} F_{i,j}}{\sum_{i \in [1, N_p]} \sum_{j \in [1, P_i]} \epsilon_{i,0} \rho_{i,0} F_{i,j}} (1 - \chi_{i,j}) \quad (3) $$

where $\epsilon_{i,0}$, $\rho_{i,0}$, and $F_{i,j}$ denote, respectively, the volume fraction of phase i, the intrinsic density of phase i and the mass fraction of sub-phase j within phase i. The subscript 0 indicates the initial time ($t = 0$).

The total production rate of species $\pi_k$ is given by

$$ \pi_k = \sum_{i \in [1, N_p]} \sum_{j \in [1, P_i]} \xi_{i,j,k} \epsilon_{i,0} \rho_{i,0} F_{i,j} \partial_t \chi_{i,j} \quad (4) $$

By summing all the contributions in the gas mixture, it is possible to obtain the overall pyrolysis-gas production rate

$$ \Pi = \sum_{k \in [1, N_g]} \pi_k \quad (5) $$

## 2.3. Mass Conservation

We assume no heterogeneous reactions between the gaseous species and the solid phase. The only production terms come from the pyrolysis.

All the solid phases, the species, and the gas mixture should be taken into account. Each of them is characterized by a different mass conservation equation. For the solid phases, the equation reads

$$ - \partial_t (\epsilon_i \rho_i) = \sum_{j \in [1, P_i]} \epsilon_{i,0} \rho_{i,0} F_{i,j} \partial_t \chi_{i,j} \quad (6) $$

For a generic species with mass fraction $y_i$, the conservation equation is

$$ \partial_t (\epsilon_g \rho_g y_i) + \partial_x \cdot (\epsilon_g \rho_g y_i v_g) + \partial_x \cdot F_i = \pi_i, \quad \forall i \in N_g^s \quad (7) $$

For the gas phase, the mass conservation needs to take into account the pyrolysis production rate. The term is added on the right hand-side

$$ - \partial_t (\epsilon_g \rho_g) + \partial_x \cdot (\epsilon_g \rho_g v_g) = - \sum_{i \in [1, N_p]} \partial_t (\epsilon_i \rho_i) = \Pi \quad (8) $$

where $v_g$ is the averaged gas velocity.

## 2.4. Momentum Conservation

The average gas velocity is obtained from the Darcy's law

$$ v_g = - \frac{1}{\epsilon_g \mu_g} \underline{K} \cdot \partial_x p_g \quad (9) $$

where $\underline{K}$ is the permeability tensor and $p_g$ the gas pressure. This expression for the velocity can be substituted back into the gas mass conservation law, Equation (8). By assuming the gas obeys the perfect gas, the overall equation becomes

$$ \partial_t \left( \frac{\epsilon_g \mathcal{M}}{RT_g} p_g \right) - \partial_x \cdot \left( \frac{p_g \mathcal{M}}{RT_g} \frac{1}{\mu_g} \underline{K} \cdot \partial_x p_g \right) = \Pi \quad (10) $$

where $\mathcal{M}$ is the mean molar mass of the gas mixture. This equation in pressure is directly solved instead of handling the mass conservation law (Equation (8)), the momentum equation (Equation (9)), and the perfect gas law.
## 2.5. Energy Conservation

The Local Thermal Equilibrium (LTE) condition is assumed between the phases: $T_s = T_g = T$. A single energy conservation is then necessary

$$ \sum_{i \in [1, N_p]} [(\epsilon_i \rho_i C_{p,i}) \partial_t T] - \partial_x \cdot (\underline{k}_{eff} \partial_x T) = - \sum_{i \in [1, N_p]} h_i \partial_t (\epsilon_i \rho_i) - \partial_t (\epsilon_g \rho_g h_g - \epsilon_g p_g) + \partial_x \cdot (\epsilon_g \rho_g h_g v_g) + \partial_x \cdot \sum_{k \in [1, N_g]} (Q_k) + \mu_g \epsilon_g^2 (\underline{K}^{-1} \cdot v_g) \cdot v_g \quad (11) $$

where $C_p$ is the heat capacity at constant pressure, $\underline{k}_{eff}$ is the effective thermal conductivity tensor, $h$ is the enthalpy, and $Q_k$ is the heat transport by effective diffusion of the species.

## 3. Numerical Model: Environment Region

The model to describe the environment region is now presented. The governing equations are derived from the Navier-Stokes equations after applying the LES filtering process and using the Favre mean variables [23]. By considering a generic $\Phi$ variable, the following notations are introduced: $\bar{\Phi}$ the filtered variable over a time period $\mathcal{T}$

$$ \bar{\Phi} = \frac{1}{\mathcal{T}} \int_{\mathcal{T}} \Phi(t) dt \quad (12) $$

and the filtered density-weighted variable $\bar{\Phi}$

$$ \bar{\Phi} = \frac{\bar{\rho}\bar{\Phi}}{\Phi} \quad (13) $$

### 3.1. Continuity Equation

The continuity equation reads

$$ \frac{\partial \bar{\rho}}{\partial t} + \nabla \cdot (\bar{\rho} \bar{\mathbf{u}}) = 0 \quad (14) $$

### 3.2. Momentum Conservation

By considering Newtonian fluids the equation takes the following form

$$ \frac{\partial \bar{\rho} \bar{\mathbf{u}}}{\partial t} + \nabla \cdot (\bar{\rho} \bar{\mathbf{u}} \bar{\mathbf{u}}) = - \nabla \bar{\rho} \bar{m} - (\mathbf{g} \cdot \mathbf{x}) \nabla \bar{\rho} + \nabla \cdot \left[ \mu_{eff} (\nabla \bar{\mathbf{u}} + (\nabla \bar{\mathbf{u}})^T - \frac{2}{3} (\nabla \cdot \bar{\mathbf{u}}) \underline{I}) \right] \quad (15) $$

where $\bar{\rho} \bar{m}$ is the modified pressure, introduced to improve the effectiveness of the numerical solution [24]. It is defined from the thermodynamic pressure as $\bar{\rho} \bar{m} = \bar{\rho} - \bar{\rho} \mathbf{g} \cdot \mathbf{x}$. The isotropic part of the sub-grid scale stress tensor is neglected since in this work Mach < 0.4 [25].

### 3.3. Energy Conservation

The equation is solved in terms of the enthalpy $\tilde{h}$

$$ \frac{\partial \bar{\rho} \tilde{h}}{\partial t} + \nabla \cdot (\bar{\rho} \tilde{\mathbf{u}} \tilde{h}) = \frac{D \bar{\rho}}{Dt} + \nabla \cdot (\bar{\rho} \alpha_{eff} \nabla \tilde{h}) + \overline{Q_{comb}} - \nabla \cdot \overline{q_R} \quad (16) $$

where $\alpha_{eff}$ is the effective thermal diffusivity, defined as the sum of the mixture thermal diffusivity $\alpha$ with the sub-grid scale thermal diffusivity $\alpha_{sgs}$. $Q_{comb}$ is the heat generated by combustion and $\overline{q_R}$ the thermal radiation flux.

### 3.4. Species Transport Equations

The gas is a mixture of different species. In order to determine its composition, a transport equation for each species is needed

$$ \frac{\partial \bar{\rho} \tilde{y}_k}{\partial t} + \nabla \cdot (\bar{\rho} \tilde{\mathbf{u}} \tilde{y}_k) = \nabla \cdot (\bar{\rho} D_{eff} \nabla \tilde{y}_k) + \overline{\omega}_i \quad (17) $$
where $D_{eff}$ is the effective mass diffusivity, that is the sum of the diffusion coefficient of each species $D_k$ with the sub-grid scale mass diffusivity $D_{sgs}$. $\omega_k$ is the rate of reaction of the generic k-th specie.

### 3.5. Ideal Gas

The gas mixture is assumed ideal. The following equations of state are considered

$$ \bar{\rho} = \frac{\bar{p}M}{R\bar{T}}, \quad M = \left[ \sum_{k=1}^{N_g} \frac{\bar{y}_k}{M_k} \right]^{-1} \qquad (18) $$

$$ \bar{h} = \sum_{k=1}^{N_g} \bar{y}_k h_k(\bar{T}), \quad h_{sk}(\bar{T}) = h_k^{ref} + \int_{T_{ref}}^{\bar{T}} C_{pk} dT \qquad (19) $$

where the Janaf model is used to evaluate the specific heat at constant pressure

$$ C_p(\bar{T}) = R \left[ \frac{a_1}{\bar{T}^2} + \frac{a_2}{\bar{T}} + a_3 + a_4 \bar{T} + a_5 \bar{T}^2 + a_6 \bar{T}^3 + a_7 \bar{T}^4 \right] \qquad (20) $$

where $a_1, \dots, a_7$ are the Janaf coefficients [26] given as input.

The dynamic viscosity $\mu$ is calculated by the Sutherland's viscosity law

$$ \mu = \frac{C_1 \bar{T}^{3/2}}{\bar{T} + C_2} \qquad (21) $$

where $C_1$ and $C_2$ are the Sutherland coefficients [27].

The thermal conductivity is determined from the modified Eucken approximation [28].

The thermal diffusivity, $\alpha$, is given by

$$ \alpha = \frac{k}{\bar{\rho}C_p}, \quad C_p = \sum_{k=1}^{N_g} \bar{y}_k C_{pk}(\bar{T}) \qquad (22) $$

### 3.6. Combustion Model

Combustion is a high-temperature exothermic chemical reaction between a fuel and an oxidant (oxygen in this case). It is modeled by means of a laminar finite rate model. The effect of turbulent fluctuations is ignored and the reaction rates are determined by Arrhenius kinetic expressions. This model is then generally inaccurate for turbulent flames and exact for laminar ones.

The generic z-th chemical reaction can be written as

$$ \sum_{k=1}^{K} \xi'_{kz} A'_{k} \rightleftharpoons \sum_{k=1}^{K} \xi''_{kz} A''_{k} \qquad (23) $$

where ' and " distinguish between the forward and backward hand-sides and K is the number of species involved in the z-th reaction. The overall reaction rate of this chemical reaction is

$$ \omega_z = k'_z \prod_k [A']_{k}^{\xi'_k} - k''_z \prod_k [A'']_{k}^{\xi''_k} \qquad (24) $$

whereas the reaction rate of the specie k is

$$ \omega_k = \frac{d[A_k'']}{dt} = \sum_z \omega_z(\xi'_k - \xi'') \qquad (25) $$

The forward and backward reaction rate constant $k'$ and $k''$ are needed in order to evaluate the reaction rates. Their calculation is done again by means of the Arrhenius law.

Finally, the heat generated by the k-th reaction is

$$ \overline{Q_{comb,z}} = \Delta h_z \overline{\omega_z} \qquad (26) $$
where $\Delta h_z$ is enthalpy variation in the z-th reaction.

### 3.7. Turbulence Model

Turbulence is modeled through the Large Eddy Simulation (LES) method. The one-equation eddy viscosity model has been used to model the sub-grid scale stress tensor. The following conservation equation for the turbulent kinetic energy is solved

$$ \frac{\partial(\bar{\rho}k)}{\partial t} + \nabla \cdot (\bar{\rho}\bar{\tilde{u}}k) = \nabla \cdot (\mu_{eff} \nabla k) - \bar{\rho} \underline{\underline{\tau}}_{sgs} : \underline{\underline{\tilde{S}}} = C_{B,1} \bar{\rho} \frac{k^{1.5}}{\Delta} \quad (27) $$

where $\Delta$ is the local filter cutoff width, calculated by taking the cubic root of the grid cell volume. $\underline{\underline{\tilde{S}}}$ is the strain rate tensor and $C_{B,1}$ is a constant which value is usually set equal to 1.048 [29]. The Boussinesq approximation is considered to model the sub-grid scale stress tensor and dynamic viscosity

$$ \underline{\underline{\tau}}_{sgs} = -2\mu_{sgs} \frac{1}{2} [\nabla \tilde{u} + (\nabla \tilde{u})^T] + \frac{2}{3} \bar{\rho} k \underline{\underline{I}} \quad (28) $$

$$ \mu_{sgs} = \rho C_{B,2} \Delta k^{0.5} \quad (29) $$

where the constant $C_{B,2}$ is usually set equal to 0.094 [29].

### 3.8. Radiation Model

In this work no radiation model has been considered. We wanted to keep the model as simple as possible in order to prove the feasibility of the numerical tool. Moreover, for the simulations considered, the fire is so small that the radiation contribution may be negligible.

## 4. Numerical Model: Interface

The coupling between the two regions is schematized in Figure 3.

<figure><img src="image_7.png" /><figcaption>Figure 3. Illustration of the numerical boundary conditions for the coupling between the two regions. Conditions for the velocity, temperature, pressure and species concentration need to be imposed.</figcaption></figure>

As can be seen, at the interface the flow updates its velocity with the one of the gases, whereas for the pressure, it is the other way around. From the thermal point of view, the coupling is done by imposing the equality of the temperatures and the normal heat fluxes at the interface. Finally, for the species concentrations, the coupling depends on the flux direction. If the flux is coming from the flow to the porous material, then the concentration of the species inside the porous domain are taken equal to the ones in the flow. In the opposite case, it is the concentration of the species in the flow which are taken equal to the ones of the gas in the porous material.
## 5. Results

Two test cases have been considered to illustrate the capability of the simulation tool. All the details are presented in the next subsections.

### 5.1. Hydrogen vs. Methane Flames

The first case considers the combustion of two gaseous species in the same domain. Two separate finite-rate combustion reactions have been considered

$$ \mathrm{CH_4} + 2\mathrm{O_2} \rightarrow \mathrm{CO_2} + 2\mathrm{H_2O} \qquad (30) $$

$$ 2\mathrm{H_2} + \mathrm{O_2} \rightarrow 2\mathrm{H_2O} \qquad (31) $$

The first reaction describe the combustion reaction for methane, CH₄, where O₂, CO₂ and H₂O indicate respectively the oxygen, carbon dioxide and the water. The second reaction is the combustion of hydrogen, H₂. The kinetics of the two reactions has been modeled with an Arrhenius type formulation. The relative input parameters are showed in Table 1.

<table><caption>**Table 1.** Kinetics parameters of the two combustion reactions. CH₄ indicates the combustion of the methane, whereas H₂ the combustion of the hydrogen.</caption><thead><tr><th></th><th>A</th><th>β</th><th>T<sub>a</sub></th></tr></thead><tbody><tr><td>CH<sub>4</sub></td><td>5.2 × 10<sup>16</sup></td><td>0</td><td>14,906</td></tr><tr><td>H<sub>2</sub></td><td>4.74 × 10<sup>12</sup></td><td>0</td><td>10,064</td></tr></tbody></table>

The test case is schematized in Figure 4.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Illustration of the test case: methane and hydrogen are introduced, at initial time, at high temperature inside a domain of air. The combustion reactions for both species are implemented.</figcaption></figure>

The domain corresponds entirely to the environment region. Air is considered as fluid. Air is initially composed of 21% of O₂ and 79% of N₂. At time t = 0, methane and hydrogen are injected in the domain at two different positions at temperature 500 K. This temperature is high enough to trigger the two combustion reactions. Two flames propagate in the domain as shown in Figure 5.

Different time intervals are considered. For each of them the flames are colored based on the temperature. The two flames are different. The one on the left, that is the one due to the hydrogen combustion, is smaller and generates more energy than the one due to the combustion of methane. This is just a confirmation of what is known from the
theory [30]. The hydrogen mixture burns faster than the one of methane and reaches higher temperatures (2500 K against 2100 K).

This simple case proves the importance of introducing different pyrolysis species in the overall problem. The flame, and so the overall simulation, strongly depends on the species taken into account.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Temperature description in the domain. Three different time steps are considered. The figures show on the left the flame generated by the hydrogen combustion and on the right the one generated by the methane combustion. The two flames are different from each other. They differ both in terms of temperature and in terms of size.</figcaption></figure>

## 5.2. Wood Log Combustion

The second case aims to simulate the combustion of a wood log. The computational domain is schematized in Figure 6.

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Numerical domain for the simulation of the burning of a wood log. A generic hardwood composition has been considered for the log. Air, with simplified composition, has been taken as flow.</figcaption></figure>

A wood log with a porosity $\epsilon_w = 0.5$ is placed at the bottom of the domain. It is modeled as being composed of 4 phases: hemicellulose, lignin, cellulose and water. The phases composition reflects a generic hardwood. Air, with the same composition as in the previous case, is surrounding it. At time $t = 0$ the bottom part of the log is heated up by a Dirichlet condition at 800 K. Due to the heat, different pyrolysis reactions are triggered. A list of them with the considered kinetic parameters can be found in Table 2.
<table><caption>**Table 2.** Pyrolysis reactions and kinetic parameters. The four species (hemicellulose, cellulose, lignin, and water) undergo different degradation processes according to the chemical reactions listed in this table.</caption><thead><tr><th>j</th><th>Hemicellulose</th><th>F<sub>1,j</sub></th><th>A<sub>1,j</sub></th><th>E<sub>1,j</sub></th><th>m<sub>1,j</sub></th><th>n<sub>1,j</sub></th></tr></thead><tbody><tr><td>1</td><td>p<sub>1,1</sub> → 0.4H<sub>2</sub> + 0.4CO + 0.1CO<sub>2</sub> + 0.05CH<sub>4</sub> + 0.05N<sub>2</sub></td><td>0.40</td><td>7.94 × 10<sup>16</sup></td><td>195,000</td><td>1</td><td>0</td></tr><tr><td>2</td><td>p<sub>1,2</sub> → 0.4H<sub>2</sub> + 0.4CO + 0.1CO<sub>2</sub> + 0.05CH<sub>4</sub> + 0.05N<sub>2</sub></td><td>0.30</td><td>1.26 × 10<sup>7</sup></td><td>106,000</td><td>1</td><td>0</td></tr><tr><td colspan="7">Cellulose</td></tr><tr><td>1</td><td>p<sub>2,1</sub> → 0.4H<sub>2</sub> + 0.4CO + 0.1CO<sub>2</sub> + 0.05CH<sub>4</sub> + 0.05N<sub>2</sub></td><td>0.75</td><td>7.94 × 10<sup>18</sup></td><td>202,650</td><td>1</td><td>0</td></tr><tr><td>2</td><td>p<sub>2,2</sub> → 0.4H<sub>2</sub> + 0.4CO + 0.1CO<sub>2</sub> + 0.05CH<sub>4</sub> + 0.05N<sub>2</sub></td><td>0.16</td><td>1.26 × 10<sup>7</sup></td><td>245,000</td><td>1</td><td>0</td></tr><tr><td colspan="7">Lignin</td></tr><tr><td>1</td><td>p<sub>3,1</sub> → 0.4H<sub>2</sub> + 0.4CO + 0.1CO<sub>2</sub> + 0.05CH<sub>4</sub> + 0.05N<sub>2</sub></td><td>0.66</td><td>6.0 × 10<sup>7</sup></td><td>120,000</td><td>1</td><td>0</td></tr><tr><td colspan="7">Water</td></tr><tr><td>1</td><td>p<sub>4,1</sub> → H<sub>2</sub>O</td><td>1</td><td>5.13 × 10<sup>6</sup></td><td>86,000</td><td>1</td><td>0</td></tr></tbody></table>

All the phases undergo different pyrolysis reactions, each characterized by different parameters. The species, once produced, percolate through the material and eventually reach the external environment. At time $t = 0.15$ s, two sparks are simulated at the two sides of the log. These sparks have the purpose of triggering the combustion. Methane and hydrogen start to burn and the fire evolves according to the concentration of the species in the air. Some snapshots of the first physical second of the simulation are collected in Figure 7.

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Starting of the fire. Six different snapshots to represent the initial evolution of the fire. At time $t = 0.15$ s two sparks are generated on the side of the wood log to trigger the combustion. A fire is generated and start propagating in the domain. By comparing the time $t = 0.50$ s with the time $t = 0.70$ s it is possible to observe the periodic behavior of the flame, the puffing effect.</figcaption></figure>

At around 30 s of physical time, the wood log is completely pyrolyzed. At around 50 s of physical time, the fire is turned off due to the lack of fuel. Results in terms of H₂ and CH₄ concentration are showed in Figure 8.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Concentration of hydrogen and methane in the simulation. It can be observed how their concentration is time dependent. At time $t = 0.1$ s the concentration of the gases in the environment region is close to zero. With time, thanks to the advancement of pyrolysis, their concentration increases. For the two gases different values of concentration are reached. At time $t = 50$ s it can be seen that the fire has consumed the two species and their concentration is again nearly zero.</figcaption></figure>

In the figure it can be appreciated how the concentrations of the gaseous species change in time and depend on the pyrolysis activity in the material region. The latter is captured by the parameter $τ$ (Equation (3)), that describes the advancement of the overall pyrolysis process. When $τ = 1$, no chemical reaction has yet occurred. When $τ = 0$ all the pyrolysis mechanisms have taken place and only char is left in the material. Some snapshots of the evolution of this parameter during the simulation are shown in Figure 9.

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Advancement of pyrolysis of the wood log during the simulation. When $τ = 1$, no chemical reactions has still occurred. When $τ = 0$ all the pyrolysis mechanisms have taken place and only char is left in the material. At time $t = 0.5$ s it can be seen how the pyrolysis process starts from the bottom of the material and at time $t = 30$ s how all the material has been pyrolysed.</figcaption></figure>

Four different time steps have been captured, and each of them focuses on the material region. It can be appreciated how the log evolves from the initial state to the final one, where only char is left. At around 30 s of physical time, all the pyrolysis reactions are completed and gas is no longer produced. Once the flame in the environment region has consumed all the fuels ($CH_4$ and $H_2$), the fire extinguishes.

A particular effect that is captured in this simulation is the so called “puffing phenomenon”. It is an instability effect of the fire due to gravity, which results in a pulsation behavior. These periodic fluctuations of the fire are repeated in time with a frequency $f$. In Figure 7 the effect can be observed by looking the time lapses $t = 0.5$ s and $t = 0.7$ s. The two snapshots both represent the end of the respective pulsation and they are quite similar. We
can deduce that for the numerical simulation, the frequency of the puffing effect is about $f_{sim} \approx 5$ Hz. Different studies attempted to study this effect and to try to correlate it with physical dimensions [31–33]. A power law has been defined to fit the data

$$f_{pl} = 1.5 \frac{1}{\sqrt{L}} \qquad (32)$$

where L in this case is the burner diameter. By considering the diameter of the fire at the base as the measure of L, it is possible to find out that this correlation gives us $f_{pl} \approx 7$ Hz. The two measures of the puffing frequency are approximately the same.

## 6. Conclusions

This study proposes a new numerical framework to investigate the combustion of pyrolysis materials. A detailed three-dimensional pyrolysis model is coupled with fireFoam, a numerical solver for combustion processes. The pyrolysis model captures the physics in the material region and describes the interaction between a multi-phase reactive material with a multi-species reactive gas mixture. The description is done at the macroscopic scale, where the governing equations are derived from upscaling theories. FireFoam is used to describe the environment region. Its governing equations are derived from the Navier-Stokes equations after applying the LES filtering process and using the Favre mean variables. The two models are coupled at the material-environment interface, where specific conditions for velocity, pressure, temperature and species concentration are imposed. Two applications are considered to illustrate the capabilities of the tool. The first one considers the combustion of two different gaseous species with the implementation of two separate finite-rate combustion reactions. Results show how the two generated flames differ between each other in terms of temperature and reaction times. The behavior of the flames is strongly influenced by the pyrolysis species considered in the model. The second case illustrates the combustion of a wood log. Results illustrate how gaseous species are produced with time through the pyrolysis of the wood log and how these species enter in the environmental region to start and sustain the flame. Thanks to the flame and to the temperature increase, the pyrolysis processes speed up and the species concentrations continuously change in time according to the rate of reactions. Once the wood is fully pyrolyzed, the fire consumes all the fuels and extinguishes by itself. The puffing effect can be observed during the evolution of the simulation: the fire exhibits periodic fluctuations that follow each other with a numerical frequency that is approximately the same as the theoretical one.

The numerical framework that was developed during this study is available in the Porous material Analysis Toolbox based on OpenFoam (PATO) released Open Source by NASA (www.pato.ac (accessed on 14 September 2021)) [34].

**Author Contributions:** Conceptualization, methodology, software, validation, formal analysis, investigation, resources, data curation, writing—original draft preparation, writing—review and editing, visualization, supervision, project administration, funding acquisition: H.S., A.A.-S., J.L., F.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by a PhD grant awarded by Arts et Métiers Institute of Technology.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All the data are included within the article.

**Acknowledgments:** The research of H.S. was sponsored by a PhD grant awarded by Arts et Métiers Institute of Technology.

**Conflicts of Interest:** The authors declare no conflict of interest.
Nomenclature

Latin Letters:

<table><tr><td><i>I</i></td><td>identity tensor</td></tr><tr><td><b>K</b></td><td>permeability tensor [m<sup>2</sup>]</td></tr><tr><td><i>S</i></td><td>strain rate tensor [s<sup>-1</sup>]</td></tr><tr><td><b>k</b></td><td>thermal conductivity tensor [W m<sup>-1</sup> K<sup>-1</sup>]</td></tr><tr><td><i>g</i></td><td>gravity field [m s<sup>-2</sup>]</td></tr><tr><td><i>x</i></td><td>position vector [m]</td></tr><tr><td><i>q</i><sub>R</sub></td><td>thermal radiation flux [J m<sup>-1</sup> s<sup>-1</sup>]</td></tr><tr><td><i>v</i></td><td>average velocity [m s<sup>-1</sup>]</td></tr><tr><td><i>a</i><sub>1</sub>, ..., <i>a</i><sub>7</sub></td><td>Janaf coefficients</td></tr><tr><td><i>A</i></td><td>chemical species</td></tr><tr><td><i>A</i></td><td>Arrhenius law pre-exponential factor</td></tr><tr><td><i>C</i><sub>1</sub>, <i>C</i><sub>2</sub></td><td>Sutherland coefficients</td></tr><tr><td><i>C</i><sub>B1</sub>, <i>C</i><sub>B2</sub></td><td>Boussinesq approximation coefficients</td></tr><tr><td><i>C</i><sub>p</sub></td><td>specific heat at constant pressure [J kg<sup>-1</sup> K<sup>-1</sup>]</td></tr><tr><td><i>D</i></td><td>mass diffusivity [m<sup>2</sup> s<sup>-1</sup>]</td></tr><tr><td><i>E</i></td><td>Arrhenius law activation energy [J mol<sup>-1</sup>]</td></tr><tr><td><i>F</i><sub>i,j</sub></td><td>mass fraction of sub-phase j in phase i</td></tr><tr><td><i>f</i></td><td>puffing frequency [Hz]</td></tr><tr><td><i>h</i></td><td>specific absolute enthalpy [J kg<sup>-1</sup>]</td></tr><tr><td><i>k</i></td><td>turbulent kinetic energy [J]</td></tr><tr><td><i>L</i>, <i>l</i></td><td>macroscopic and microscopic characteristic lengths [m]</td></tr><tr><td><i>M</i></td><td>mean molar mass [kg mol<sup>-1</sup>]</td></tr><tr><td><i>m</i>, <i>n</i></td><td>Arrhenius law parameters</td></tr><tr><td><i>N</i><sub>g</sub>, <i>N</i><sub>p</sub></td><td>number of gaseous and solid species</td></tr><tr><td><i>p</i></td><td>average pressure [Pa]</td></tr><tr><td><i>P</i><sub>i</sub></td><td>number of sub-phases in the solid phase i</td></tr><tr><td><i>Q</i></td><td>heat flux [J m<sup>-2</sup> s<sup>-1</sup>]</td></tr><tr><td><i>R</i></td><td>perfect gas constant [J mol<sup>-1</sup> K<sup>-1</sup>]</td></tr><tr><td><i>s</i></td><td>generic solid phase</td></tr><tr><td><i>T</i></td><td>average temperature [K]</td></tr><tr><td><i>T</i></td><td>time period [s]</td></tr><tr><td><i>y</i></td><td>species mass fraction</td></tr></table>

Greek Letters:

<table><tr><td><i>&tau;</i></td><td>stress tensor [N m<sup>-2</sup>]</td></tr><tr><td><i>&alpha;</i></td><td>thermal diffusivity [m<sup>2</sup> s<sup>-1</sup>]</td></tr><tr><td><i>&epsilon;</i></td><td>volume fraction</td></tr><tr><td><i>&mu;</i></td><td>dynamic viscosity [Pa s]</td></tr><tr><td><i>&zeta;</i></td><td>mass stoichiometric coefficient</td></tr><tr><td><i>&pi;</i><sub>k</sub></td><td>pyrolysis production rate of species k [kg m<sup>-3</sup> s<sup>-1</sup>]</td></tr><tr><td><i>&Pi;</i></td><td>total pyrolysis gas production rate [kg m<sup>-3</sup> s<sup>-1</sup>]</td></tr><tr><td><i>&rho;</i></td><td>average density [kg m<sup>-3</sup>]</td></tr><tr><td><i>&tau;</i></td><td>overall pyrolysis advancement</td></tr><tr><td><i>&chi;</i><sub>1,j</sub></td><td>advancement of pyrolysis reaction j within phase i</td></tr><tr><td><i>&omega;</i><sub>k</sub></td><td>combustion rate of reaction of the species k [kg m<sup>-3</sup> s<sup>-1</sup>]</td></tr></table>

Subscripts and Superscripts:

<table><tr><td>comb</td><td>combustion</td></tr><tr><td>eff</td><td>effective</td></tr><tr><td><i>f</i></td><td>fluid</td></tr><tr><td><i>g</i></td><td>gas</td></tr><tr><td><i>i</i></td><td>index for the solid phases</td></tr><tr><td><i>j</i></td><td>index for the sub-phases produced from pyrolysis</td></tr><tr><td><i>k</i></td><td>index for the gaseous species</td></tr></table>
<table><tr><td>m</td><td>modified</td></tr><tr><td>pl</td><td>power law</td></tr><tr><td>s</td><td>solid</td></tr><tr><td>sgs</td><td>sub-grid scale</td></tr><tr><td>sim</td><td>simulation</td></tr><tr><td>w</td><td>wood</td></tr><tr><td>z</td><td>index for the combustion chemical reactions</td></tr><tr><td colspan="2"><strong>Acronyms:</strong></td></tr><tr><td>LES</td><td>Large Eddy Simulation</td></tr><tr><td>LTE</td><td>Local Thermal Equilibrium</td></tr><tr><td>PATO</td><td>Porous material Analysis Toolbox based on OpenFoam</td></tr><tr><td>REV</td><td>Representative Elementary Volume</td></tr></table>

## References

1. Östman, B.; Brandon, D.; Frantzich, H. Fire safety engineering in timber buildings. *Fire Saf. J.* **2017**, *91*, 11–20. [CrossRef]

2. Cox, G.; Moss, J. Fire science and aircraft safety. In Proceedings of the AGARD Conference Proceedings, The Hague, The Netherlands, 8–12 May 1989.

3. Mouritz, A.P.; Gibson, A.G. *Fire Properties of Polymer Composite Materials*; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2007; Volume 143.

4. Dasari, A.; Yu, Z.Z.; Cai, G.P.; Mai, Y.W. Recent developments in the fire retardancy of polymeric materials. *Prog. Polym. Sci.* **2013**, *38*, 1357–1387. [CrossRef]

5. Delichatsios, M.; Paroz, B.; Bhargava, A. Flammability properties for charring materials. *Fire Saf. J.* **2003**, *38*, 219–228. [CrossRef]

6. Fateh, T.; Richard, F.; Batiot, B.; Rogaume, T.; Luche, J.; Zaida, J. Characterization of the burning behavior and gaseous emissions of pine needles in a cone calorimeter—FTIR apparatus. *Fire Saf. J.* **2016**, *82*, 91–100. [CrossRef]

7. White, R.H.; Dietenberger, M.A. Fire safety. In *Wood Handbook: Wood as an Engineering Material*; General Technical Report FPL; USDA Forest Service, Forest Products Laboratory: Madison, WI, USA, 1999; Volume 113, pp. 17.1–17.16.

8. Global, F. fireFoam. Available online: http://code.google.com/p/firefoam-dev/ (accessed on 14 September 2021).

9. Li, Y.Z.; Huang, C.; Anderson, J.; Svensson, R.; Ingason, H.; Husted, B.; Runefors, M.; Wahlqvist, J. Verification, Validation and Evaluation of fireFoam as a Tool for Performance Design. 2017. Available online: http://www.diva-portal.org/smash/record.jsf?pid=diva2%3A1166638 (accessed on 14 September 2021).

10. Zamorano, R. fireFoam (CFD Solver) Validation in Compartment Fire Scenario Using High Resolution Data. 2018. Available online: https://wwwsemanticscholar.org/paper/FireFoam-(CFD-solver)-validation-in-compartment-Zamorano/bd0ba6c49d5073adc6379008b56ec9182c1ec164 (accessed on 14 September 2021).

11. Trouvé, A.; Wang, Y. Large eddy simulation of compartment fires. *Int. J. Comput. Fluid Dyn.* **2010**, *24*, 449–466. [CrossRef]

12. Yang, H.; Yan, R.; Chen, H.; Lee, D.; Zheng, C. Characteristics of hemicellulose, cellulose and lignin pyrolysis. *Fuel* **2007**, *86*, 1781–1788. [CrossRef]

13. Blondeau, J.; Jeanmart, H. Biomass pyrolysis at high temperatures: Prediction of gaseous species yields from an anisotropic particle. *Biomass Bioenergy* **2012**, *41*, 107–121. [CrossRef]

14. Nunn, T.; Howard, J.; Longwell, J.; Peters, W. Product compositions and kinetics in the rapid pyrolysis of milled wood lignin. *Ind. Eng. Chem. Process. Des. Dev.* **1985**, *24*, 844–852. [CrossRef]

15. Di Blasi, C. Modeling chemical and physical processes of wood and biomass pyrolysis. *Prog. Energy Combust. Sci.* **2008**, *34*, 47–90. [CrossRef]

16. Chan, W.; Kelbon, M.; Krieger, B. Modelling and experimental verification of physical and chemical processes during pyrolysis of a large biomass particle. *Fuel* **1985**, *64*, 1505–1513. [CrossRef]

17. Miller, R.; Bellan, J. A generalized biomass pyrolysis model based on superimposed cellulose, hemicellulose and lignin kinetics. *Combust. Sci. Technol.* **1997**, *126*, 97–137. [CrossRef]

18. Ranzi, E.; Cuoci, A.; Faravelli, T.; Frassoldati, A.; Migliavacca, G.; Pierucci, S.; Sommariva, S. Chemical kinetics of biomass pyrolysis. *Energy Fuels* **2008**, *22*, 4292–4300. [CrossRef]

19. Lachaud, J.; Magin, T.; Cozmuta, I.; Mansour, N. A short review of ablative-material response models and simulation tools. In Proceedings of the 7th European Symposium on Aerothermodynamics, SP-692, European Space Agency, Noordwijk, The Netherlands, 9–12 May 2011; pp. 1–12.

20. Lachaud, J.; Mansour, N. Porous-material analysis toolbox based on OpenFOAM and applications. *J. Thermophys. Heat Transf.* **2014**, *28*, 191–202. [CrossRef]

21. Lachaud, J.; van Eekelen, T.; Scoggins, J.; Magin, T.; Mansour, N. Detailed chemical equilibrium model for porous ablative materials. *Int. J. Heat Mass Transf.* **2015**, *90*, 1034–1045. [CrossRef]

22. Lachaud, J.; Scoggins, J.; Magin, T.; Meyer, M.; Mansour, N. A generic local thermal equilibrium model for porous reactive materials submitted to high temperatures. *Int. J. Heat Mass Transf.* **2017**, *108*, 1406–1417. [CrossRef]

23. Favre, A. Turbulence: Space-time statistical properties and behavior in supersonic flows. *Phys. Fluids* **1983**, *26*, 2851–2863. [CrossRef]
24. Ferziger, J.H.; Perić, M.; Street, R.L. *Computational Methods for Fluid Dynamics*; Springer: Berlin/Heidelberg, Germany, 2002; Volume 3.

25. Erlebacher, G.; Hussaini, M.Y.; Speziale, C.G.; Zang, T.A. Toward the large-eddy simulation of compressible turbulent flows. *J. Fluid Mech.* **1992**, *238*, 155–185. [CrossRef]

26. Chase, M. NIST-JANAF Thermochemical Tables. 1998. Available online: https://janaf.nist.gov/ (accessed on 14 September 2021).

27. Shen, V.K.; Siderius, D.W.; Krekelberg, W.P.; Hatch, H.W. NIST Standard Reference Simulation Website-SRD 173. *NIST Standard Reference Database Number 173*; National Institute of Standards and Technology: Gaithersburg, MD, USA, 2017.

28. Poling, B.E.; Prausnitz, J.M.; O'Connell, J.P. *Properties of Gases and Liquids*; McGraw-Hill Education: New York City, NY, USA, 2001.

29. Xiao, H.; Jenny, P. A consistent dual-mesh framework for hybrid LES/RANS modeling. *J. Comput. Phys.* **2012**, *231*, 1848–1865. [CrossRef]

30. Sierens, R.; Demuynck, J.; Paepe, M.D.; Verhelst, S. Heat Transfer Comparison between Methane and Hydrogen in a Spark Ignited Engine. 2010. Available online: https://juser.fz-juelich.de/record/135724/files/TA3_3_Sierens.pdf (accessed on 14 September 2021).

31. Mandin, P.; Most, J. Characterization of the puffing phenomenon on a pool fire. *Fire Saf. Sci.* **2000**, *6*, 1137–1148. [CrossRef]

32. Beyler, C.L. Fire plumes and ceiling jets. *Fire Saf. J.* **1986**, *11*, 53–75. [CrossRef]

33. Pretrel, H.; Audouin, L. Periodic puffing instabilities of buoyant large-scale pool fires in a confined compartment. *J. Fire Sci.* **2013**, *31*, 197–210. [CrossRef]

34. Lachaud, J.; Meurisse, J. Porous Material Analysis Toolbox Based on OpenFOAM. Available online: http://pato.ac/ (accessed on 14 September 2021).