Simulation of Wood Combustion in PATO Using a Detailed Pyrolysis Model Coupled to ﬁreFoam

: The numerical simulation of ﬁre 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 ﬂame, which produces the energy that sustains the pyrolysis process. Experimental and numerical works of the ﬁre 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 ﬁre propagation, we have developed a numerical framework that couples a detailed three-dimensional pyrolysis model and ﬁreFoam. 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 ﬁrst proof of concept of the coupling of a multi-species pyrolysis process and a ﬂame. In the ﬂame, we consider two separate ﬁnite-rate combustion reactions for methane and hydrogen. The ﬂame 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 ﬂame.


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][10][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 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.  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][21][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.

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.

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.

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 ξ i,j,k as follows 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 χ i,j of sub-phase j within phase i can be defined as where R is the perfect gas constant, T g the temperature of the gas, A the Arrhenius law pre-exponential factor, E the Arrhenius law activation energy, and m and n the Arrhenius law parameters. The total advancement in the pyrolysis process, τ, is evaluated as follows where i , ρ i , 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 π k is given by By summing all the contributions in the gas mixture, it is possible to obtain the overall pyrolysis-gas production rate

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 For a generic specie with mass fraction y i , the conservation equation is 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 where v g is the averaged gas velocity.

Momentum Conservation
The average gas velocity is obtained from the Darcy's law where 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 where 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.

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 where C p is the heat capacity at constant pressure, k e f f is the effective thermal conductivity tensor, h is the enthalpy, and Q k is the heat transport by effective diffusion of the species.

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 Φ variable, the following notations are introduced: Φ the filtered variable over a time period T and the filtered density-weighted variableΦ

Continuity Equation
The continuity equation reads

Momentum Conservation
By considering Newtonian fluids the equation takes the following form where p m is the modified pressure, introduced to improve the effectiveness of the numerical solution [24]. It is defined from the thermodynamic pressure as p m = p − ρ g · x.
The isotropic part of the sub-grid scale stress tensor is neglected since in this work Mach < 0.4 [25].

Energy Conservation
The equation is solved in terms of the enthalpyh ∂ρh ∂t where α e f f is the effective thermal diffusivity, defined as the sum of the mixture thermal diffusivity α with the sub-grid scale thermal diffusivity α sgs . Q comb is the heat generated by combustion and q R the thermal radiation flux.

Species Transport Equations
The gas is a mixture of different species. In order to determine its composition, a transport equation for each specie is needed where D e f f 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 . ω k is the rate of reaction of the generic k-th specie.

Ideal Gas
The gas mixture is assumed ideal. The following equations of state are considered where the Janaf model is used to evaluate the specific heat at constant pressure where a 1 , .., a 7 are the Janaf coefficients [26] given as input.
The dynamic viscosity µ is calculated by the Sutherland's viscosity law 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, α, is given by

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 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 whereas the reaction rate of the specie k is 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 where ∆h z is enthalpy variation in the z-th reaction.

Turbulence Model
Turbulence is modeled through the Large Eddy Simulation (LES) method. The oneequation 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 where ∆ is the local filter cutoff width, calculated by taking the cubic root of the grid cell volume.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 µ sgs = ρC B,2 ∆k 0.5 (29) where the constant C B,2 is usually set equal to 0.094 [29].

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.

Numerical Model: Interface
The coupling between the two regions is schematized in Figure 3. 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.

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

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 The first reaction describe the combustion reaction for methane, CH 4 , where O 2 , CO 2 and H 2 O indicate respectively the oxygen, carbon dioxide and the water. The second reaction is the combustion of hydrogen, H 2 . The kinetics of the two reactions has been modeled with an Arrhenius type formulation. The relative input parameters are showed in Table 1. The test case is schematized in Figure 4. The domain corresponds entirely to the environment region. Air is considered as fluid. Air is initially composed of 21% of O 2 and 79% of N 2 . 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.

Wood Log Combustion
The second case aims to simulate the combustion of a wood log. The computational domain is schematized in Figure 6.  A wood log with a porosity 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. 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. 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 2 and CH 4 concentration are showed in Figure 8. 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. 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 laps 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 ≈ 5 Hz. Different studies attempted to study this effect and to try to correlate it with physical dimensions [31][32][33]. A power law has been defined to fit the data 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 ≈ 7 Hz. The two measures of the puffing frequency are approximately the same.

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].