Within the gaseous domain, the Favre-filtered compressible Navier-Stokes equations are solved along with conservation equations for mass, total enthalpy, and reacting species. These equations are given as (Lapointe et al., 2021)

where ρ is density, U is velocity, p_rgh is dynamic pressure, p is total pressure, g is gravitational acceleration, h is specific sensible enthalpy, K is specific kinetic energy, and Y_i is the mass fraction of the i^th specie. Unity Lewis and Prandtl numbers are assumed and a one equation eddy-viscosity turbulence model is used to represent the turbulent contributions to the effective viscosity and thermal diffusivity, denoted μ_eff and α_eff, respectively. The heat transfer terms Q_rxn and Q_rad represent effects due to reactions and radiations, respectively, and ω_i is the reaction rate of the ith specie.

Within the solid, one-dimensional conservation equations for mass, reacting species, and enthalpy are solved to model heat transfer and pyrolysis (Chaos et al., 2011; Vinayak, 2017; Fukumoto et al., 2018). These equations are given as (Lapointe et al., 2021)

where ρ_s is the solid density, Ṙ_gas is pyrolysis gas production rate, Y_k,s is the mass fraction of the k^th solid chemical specie, Ṙ_k is the corresponding reaction rate, h_s is the solid-phase specific enthalpy, and κ_s is the thermal conductivity. Source terms are included in the solid-phase enthalpy equations to account for radiation (S_rad), reactions (S_rxn), gas production (S_gas), and gas motion within the solid (S_flux).

Gas-phase combustion is modeled using infinitely-fast chemistry. Prior work with fireDyMFoam has shown that infinitely-fast chemistry can be used to accurately reproduce experimental measurements of diffusion flames in pool fires (Lapointe et al., 2020) and during solid-phase pyrolysis (Lapointe et al., 2021), supporting the use of this computational simplification in the present study.

The solid- and gas-phase regions are coupled using mixed, mapped boundary conditions for temperature, fuel mass fraction, and velocity. Temperature and velocity boundary conditions follow those described by Vinayak (2017). Fuel mass fractions are modeled using a look-up table approach to assign temperature-dependent compositions to the gas-phase pyrolysis products, as described in the next section.

From a physical standpoint, the pyrolysis process consists of thousands to millions of elementary reactions (Demirbaş, 2000). Reduced-order models are thus required to render the pyrolysis process computationally tractable. Attempts to capture the complexity of the pyrolysis mechanism with significantly simpler kinetics have included single-component/single-step, single-component/multi-step, multi-component/single-step, and multi-component/multi-step kinetics (Di Blasi, 2008). Ranzi et al. (2008) published a gas-phase chemical mechanism that relied on species that are produced in their multi-component/multi-step pyrolysis kinetic model. These models have been updated over the years with the assistance of pyrolysis experiments, notably by Corbetta et al. (2014), Debiagi et al. (2015), and Gentile et al. (2017). One common factor of these works is the lack of gaseous ignition or combustion modeling. The present study aims to combine the existing fireFoam (and by extension fireDyMFoam) framework for modeling solid-phase pyrolysis and gas-phase combustion with the multi-component/multi-step pyrolysis model published by Debiagi et al. (2015).

For simplicity and computational speed, we predict the rate of formation of char and gas-phase pyrolysis products using a single-component, single-step approach with an irreversible Arrhenius reaction, as given by

with the pyrolysis reaction rate, k_pyro, given by

Here, n is the order of the reaction, A is the pre-exponential factor, β is the temperature dependent constant, T is the local solid-phase temperature, T_a is the activation temperature, where T_a ≡ E_a/R_u (namely, the activation energy divided by the universal gas constant). Temperature-varying specific heats are applied for wood and char in the form of a power-law fit given by

where c₀, n₀ and T_ref are reference values of the specific heat, exponent, and temperature that must be defined for each fuel sample.

As defined in Eq. (9), the solid-phase reaction only proceeds after the critical temperature T_crit is reached. Once produced, the pyrolysis gases are immediately transferred out of the solid; this approximation was recently validated by Agarwal et al. (2021). In the present study, we assign a temperature-dependent composition of these gases using a custom boundary coupling approach to convert the generic “gas” produced by pyrolysis of the solid-phase to multi-species products predicted by Debiagi et al. (2015). This study uses multi-step, multi-species kinetics to convert the main constituents of wood (i.e., cellulose, hemicellulose, lignin, and extractives) to combustible materials, intermediates, and char in a zero-dimensional process.

To incorporate this model within the OpenFOAM-7 pyrolysis framework, the equilibrium species mass fractions are tabulated as a function of temperature from 300 to 1,500 K in 5 K increments. Calculations were performed zero-dimensionally (i.e., as functions of time only) using ordinary differential equation integration (ode23s) in MATLAB, and the composition corresponding to each temperature was determined from the steady-state, or equilibrium, mass fractions. These are applied as look-up table functions at the solid-gas boundary in a custom version of the totalFlowRateAdvectiveDiffusive boundary condition. Here, the average temperature across the solid-gas interface is used to evaluate the pyrolysate composition by relative mass fractions of species entering the gaseous domain. Examples of these boundary conditions are included in the sample case found in the public github repository where fireDyMFoam can also be found (github.com/clapointe2011/public). While these new boundary conditions rely on the average surface temperature of the solid, they could also be extended to a cell-by-cell basis depending on the cell-valued temperature.

In the gas domain, first-order temporal integration and second-order spatial discretization are used. Total variation diminishing variants of the central differencing scheme are used for scalar divergence terms, stabilized central differencing is used for velocity divergence, and enthalpy gradients are limited to bound temperature. The solid region is numerically configured following cases provided with OpenFOAM (The OpenFOAM Foundation, 2019), resulting in first-order in time and second-order in space accuracy. Maximum advective Courant and solid diffusion numbers of 0.4 and 1, respectively, are used to set a global time step. Time integration is performed using the pressure-implicit with splitting of operators (PISO) algorithm, and a fixed number of 3 PISO loops are used. We use a finite-volume implementation of the discrete ordinate method to model radiative heat transfer, with the grey mean absorption emission radiation model, prescribed coefficients based on fireFoam tutorials, and 100 discrete angles.

Load-balanced AMR is used in the gas-phase region, where refinement is based on multiple fields simultaneously. Further details on the specific configuration of the AMR used here is provided in the next section, and additional description of the development, verification, and validation of AMR in the fireDyMFoam solver is provided in Lapointe et al. (2020, 2021).

The simulations use a 1 m³ total domain size and a [0.12 × 0.04 × 0.02] m Douglas fir solid fuel sample to represent the experiments. Although this simulation configuration, shown in Figure 2, is an idealization of the experimental configuration provided schematically in Figure 1, the available experimental measurements were all taken at or near the surface of the Douglas fir sample, away from the cone-calorimeter, structural supports, and other diagnostic equipment. Thus, the influence of the surroundings is minimal and it is not necessary to model the full experimental geometry to obtain good agreement with the measurements, as will be shown in Section 4.4.

The gaseous domain is initially filled with a quiescent mixture of air and water vapor to match the ambient experimental humidity. All sides of the gaseous region are open and allow for entrainment, with the exception of the base, which is modeled as a mild air co-flow (5 mm/s), and the solid-gas interface. This coupling takes the form of boundary conditions for velocity (to conserve the mass flux of gas from the solid region to pyrolysate species in the gaseous region), temperature (to conserve energy across the gas-solid interface), and species (converting generic “gas” product in the solid region to known species in the gaseous region). At 89 s, combustion is initiated via infinitely-fast chemistry, matching the time when the experiment was ignited via a spark.

In the solid region, physical properties are set to match those of the Douglas fir sample used in the experiment. Properties of the wood are assumed to be homogeneous, with wood density ρ_wood = 524 kg/m³ based on experimental measurements, char density ρ_char = 73 kg/m³, thermal conductivities κ_wood = 0.11 W/(m·K) and κ_char = 0.065 W/(m·K), emissivities of 0.759 and 0.957 for raw and charred wood, respectively. where c₀, n₀ and T_ref are given in Table 1, and the specific heat of formation is set as hs,wood=-4×104 kJ/kg (Parker, 1989; Ragland et al., 1991). To model the radiative heat source used in the experiment, a heat flux of 21.9 kW/m² is applied uniformly to the interface following the implementation in Vinayak (2017). The variable specific heat of Douglas fir wood and char are prescribed by parameters in Table 1. The pyrolysis reaction parameters are those in Table 2.

Initial mole fractions of cellulose, hemicellulose, LIG-C, LIG-H, LIG-O, TANN, and TGL for Douglas Fir from Debiagi et al. (2015), Faravelli et al. (2010), and Schwetz and Lipp (1985); Table 3 lists these in terms of weight percentages. Previous work (Glusman et al., 2019) used the experimental surface temperature measurements of Makowiecki et al. (2020a) to calculate the pyrolysate composition and compare to measured water mole fraction.

The composition of the gas-phase pyrolysis products corresponding to the Douglas fir sample are obtained by running out isothermal pyrolysis to steady state at various temperatures, as described in Section 3.2. The resulting concentrations of the eighteen pyrolysate species are shown as functions of temperature in Figure 3, and are constrained such that the sum over all species mass fractions is unity.

The computational domain is discretized at the coarsest level with a 4-cm uniform mesh. The cells within a [0.20 × 0.25 × 0.10] m region centered on the bottom boundary are refined by three levels (i.e., eight times finer in each coordinate direction), corresponding to the region of interest. The solid-gas interface containing cells within a [0.12 × 0.04] m rectangle, centered at the bottom of the domain, are refined with an additional level of AMR and are extruded downwards by 0.02 m to create the computational mesh for the solid. The resulting gas region mesh is composed of roughly 40,000 cells. The solid-phase domain has 192 one-dimensional (vertically oriented) regions, each with 45 cells, for 8,640 total cells.

Within the gaseous region, three levels of AMR are used to attain 5 mm resolution in areas of active mixing of air and pyrolysate. Refinement is performed in areas of high methane concentration and high gradients of methane concentration. Heat release is also included in the AMR criteria for refining areas of active combustion. Each parameter is scaled between values of zero and one; using methane concentration as an example, the normalized concentration is given by

where [CH4^] is the normalized methane concentration, and the absolute methane concentration is [CH₄]. The criteria for refinement is based on the top 97.5% of normalized methane concentration, the top 99% of normalized methane concentration gradient, and the top 99.9% of normalized heat release.

An example AMR mesh at one time step, including the corresponding temperature field, is shown in Figure 4. Although the mesh does not change substantially during the later stages of the simulation, the AMR effectively provides a minimal mesh required to maintain the target resolution, resulting in an efficient numerical approach that does not require a priori knowledge of the flow. The majority of the simulation is run with roughly 65,000 cells. Comparable results were achieved with static mesh refinement (SMR) with nearly 385,000 cells.

Because data collected from the frequency combs is line-of-sight averaged, simulated values for gas-phase temperature and species must be similarly processed. Values are linearly interpolated onto horizontal profiles across the center of the block at 2, 4, and 6 mm above the surface due to uncertainty in laser location attributable to swelling of the wood surface and beam steering (due to changes in the index of refraction) during the course of the experiment. Profile averages are then compared to experimental data.

Good agreement is observed between experimental and computational results for the solid surface temperature, as shown in Figure 5 for the simulated cases. The experimental surface temperatures were determined from the infrared camera using a constant surface emissivity, a time-varying mixture-dependent (wood and char) value is used computationally. The experimental mean was calculated from the beginning of each run until one experiment ignited, after that a dotted line is used to indicate an extrapolated experimental mean. Post-ignition surface temperature measurements are not shown for the experiment as the flame interferes with the measurement accuracy. Future experiments will utilize embedded thermocouples to capture the surface temperature post-ignition at near-surface locations and subsequent depths to further verify the solid model.

Figure 6 indicates that the simulated cumulative mass loss is accurate in capturing the mass loss during most of the combustion process. The experimental and simulated data are centered with respect to one another by setting the time of ignition to t = 0 s. Potential sources of discrepancy are that char oxidation and the physics of smoldering combustion are not explicitly solved in this model. Furthermore, the one-dimensional nature of the solid region inhibits heat transfer transversely and thus diminishes the mass loss at later times where high surface temperatures are prevalent. Of particular interest, even in the absence of finite-rate chemistry in the gas-phase combustion, it is possible to obtain agreement with the solid-surface temperature and cumulative mass loss with the time-resolved experiments. The mass loss rates from the simulations and experiments are shown in Figure 7, revealing reasonable agreement during the pyrolysis phase and qualitative agreement of the mass loss trends post-ignition. These data are also ignition-centered in time.

The simulated time-series of gas-phase temperature is shown in Figure 8. Included in the figure are “clouds” of simulated temperatures at the nominal height of 4 ± 2 mm to capture uncertainty in the measurement height. Two factors cause the height uncertainty: as pyrolysis occurs, it is expected that the wood sample surface would swell and deform in-homogeneously and laser beam steering due to the index of refraction changing in the gas above the sample. The gray lines in Figure 8 show results from seven independent experiments of the same configuration. These results are centered at the time of ignition, enabling easier identification of temporal variations pre- and post-ignition. Simulations and averaging were carried out with the average ignition time of 89 s. Figure 8 shows that the simulation produces good agreement of gas-phase temperature after ignition although the gas-phase temperature is over-predicted during pyrolysis. This indicates that more sophisticated modeling of the solid-phase may be necessary, which is an important direction for future research.

The water vapor mole fraction, normalized by ignition time, is shown in Figure 9 along with “clouds” of simulated temperatures at the nominal height of 4 ± 2 mm to capture uncertainty in the measurement height. A deficit of water vapor during pyrolysis and the deficit persists after ignition. This discrepancy points to the need for additional improvements to the pyrolysis boundary conditions post-ignition and to the composition of the pyrolysate. However, given the simplicity (i.e., infinitely-fast chemistry and use of only methane and water as the pyrolysate) of the computational model used here, the present results are promising for future simulations.