A low-Mach, unstructured, large-eddy-simulation-based, unsteady flamelet approach with a generalized heat loss combustion methodology (including soot generation and consumption mechanisms) is deployed to support a large-scale, quiescent, 5-m JP-8 pool fire validation study. The quiescent pool fire validation study deploys solution sensitivity procedures, i.e., the effect of mesh and time step refinement on capturing key fire dynamics such as fingering and puffing, as mesh resolutions approach O(1) cm. A novel design-order, discrete-ordinate-method discretization methodology is established by use of an analytical thermal/participating media radiation solution on both low-order hexahedral and tetrahedral mesh topologies in addition to quadratic hexahedral elements. The coupling between heat losses and the flamelet thermochemical state is achieved by augmenting the unsteady flamelet equation set with a heat loss source term. Soot and radiation source terms are determined using flamelet approaches for the full range of heat losses experienced in fire applications including radiative extinction. The proposed modeling and simulation paradigm are validated using pool surface radiative heat flux, maximum centerline temperature location, and puffing frequency data, all of which are predicted within 10% accuracy. Simulations demonstrate that under-resolved meshes predict an overly conservative radiative heat flux magnitude with improved comparisons as compared to a previously deployed hybrid Reynolds-averaged Navier–Stokes/eddy dissipation concept-based methodology.

## I. INTRODUCTION

Large-scale pool fires are of interest to applications that involve the transport or storage of sensitive resources that are in close proximity to hydrocarbon-based liquid fuels. Accident scenarios involving large inventory loss of liquid fuel support the existence of pool fires. Peak convective and radiative loads are a function of pool fire size, shape, and crosswind and the damaged state of the transportation device and surroundings.^{1,2} In most cases, detailed understanding of heat fluxes to both near and far objects (both large and small) is classically of interest.^{3} In the most general description, a pool fire is represented by a tightly coupled physics set including buoyant turbulent plumes, gas-phase reactions, soot generation and consumption, and participating media radiation (PMR).^{4} In most cases, conjugate heat transfer effects are required. Therefore, radiative fluxes in and around the pool fire are of general interest to both predicting devolatilization of fuel and quantifying thermal insults to valuable objects that are either fully engulfed or in proximity to the fire.^{5}

Sandia National Laboratories holds a rich history of developing both experimental techniques and datasets, with deployment of computational approaches that are required to model pool fires.^{4} The first-generation of pool fire simulation approaches at Sandia rested on hybrid-Reynolds-Averaged Navier–Stokes (RANS) approaches in the context of Eddy Dissipation Concept (EDC) combustion modeling with empirically based soot and soot-precursor models.^{4,6} This RANS-based strategy mimicked the turbulence and combustion methods found in many early fire modeling approaches.^{7–9} Several model validation studies for this first-generation model suite exist, including non-reacting helium plumes using a hybrid-RANS approach,^{10} hydrogen- and methane-based fires,^{11} and scenarios including crossflow.^{12,13} To showcase the effort toward extreme-scale pool fire experimental characterization, pool fire experiments approaching 100 m have been conducted.^{14}

Primary interest in pool fire modeling is the desire to capture intermittent fire dynamics features such as puffing and fingering.^{5,15} Moreover, the complexity of the fire environment can substantially increase via the presence of crosswind that drives large-scale column vortex formation and increased mixing that supports more efficient combustion efficiency (especially with obstructions).^{1,16} Each of these requirements on turbulence models has realized a required transition from a RANS-based approach to a large-eddy simulation (LES) approach.^{17,18} Modern fire-based simulation tools almost exclusively employ the LES technique (see the National Institute of Standards Fire Dynamics Simulator, FDS^{19} or the Open-FOAM-based tool, Fire-FOAM^{20,21}) while in most cases retaining the EDC combustion model. More recently, flamelet-based combustion models have been introduced with sufficient fidelity to predict a range of fire-relevant phenomena.^{22–24}

To support a transition from a RANS/EDC- to a LES/flamelet-based high-fidelity approach, recent efforts within our research group have established design-order, low-dissipation numerics suitable for complex configurations that support hybrid meshing procedures.^{25} Several canonical validation studies^{25} that include structural uncertainty quantification^{26,27} have been deployed for the underlying design-order numerical methodology.^{28} Each of these works has been designed to support the fire use case while establishing capability-class, i.e., greater than 1000 nodes, simulation tools.^{29,30}

Although LES-based pool fire simulation method development and corresponding experimental studies have been conducted in the small regime, i.e., less than 1 m, few studies exist for the large-scale accident scenarios of interest to this project. For example, in the study by Weckman and Strong,^{31} a medium-scale (31-cm) methanol pool fire experiment was used to drive foundational understanding of the complex three dimensional pool fire dynamics, i.e., existence of both puffing and vortical structures. In the works of Tieszen *et al.*,^{4} the complex puffing phenomena based on arguments from the vorticity transport equation were provided and later reviewed by Tieszen.^{5} The RANS-based vorticity generation modeling approach of Nicolette *et al.*^{32} displayed limited success due to the fact that with mesh refinement, a procedure designed to capture the effect of numerical errors, transient behavior was predicted (a violation of the ergodic averaging hypothesis).

In a recent LES validation study for the Weckman and Strong pool fire configuration,^{31} a configuration that does not require a soot model, the EDC-based approach for this medium-scale methanol pool fire was conducted using FireFOAM,^{20} while several LES-based approaches have studied the effect of the experimental pool configuration.^{33,34} In the study by Koo *et al.*,^{35} soot model improvements were outlined within the LES mixture fraction-based paradigm with focus on the end pool fire use case. In the study by Yao *et al.*,^{36} ethylene/propane pool fire modeling was presented using a tunable soot onset and oxidation constants that inform Heaviside functions on soot rates on small-scale pool diameters (10-cm). Large-scale fire interaction effects that include aircraft crash-and-burn scenarios using LES were captured by Wang and Wang,^{18} who again incorporated an EDC-based combustion model. Other extensions to EDC-based combustion can be found in the study by Chen *et al.*^{37} A high-fidelity LES-based methodology was established by the Center for Accidental Fires and Explosions,^{38} where a structured, second-order in space-and-time discretization methodology was used in concert with a Gibbs free-energy minimization combustion model.^{39} Accident scenarios including pool diameters that approach the engineering scale of interest, specifically greater than 5 m, and up to 10 m were conducted as part of the Center's modeling campaign.

Although some of the smaller scale experiments (both reacting and non-reacting) conducted at Sandia included detailed experimental diagnostics such as particle image velocimetry^{11} or planar laser-induced fluorescence,^{40} the ability to measure velocity, temperature, soot, and other quantities of interest at the 5 m and beyond scale remains a significant challenge. As such, limited experimental data points for pool radiative heat flux and temperatures^{2,15} exist at the large-scale size of interest. Consequently, experimentally validated correlations for the plume height, which can be found in the study by Heskestad,^{41} are of interest to exercise within a mesh refinement modeling and simulation validation construct.

In addition to limited experimental data sets in the large-scale regime, several other complexities exist, most notably the computational scale required to adequately capture the most rudimentary pool fire observations such as core collapse. The three primary mechanisms that drive overall quiescent pool fire dynamics include (1) mis-alignment of density and hydrostatic pressure gradients, i.e., baroclinic torque, (2) Rayleigh/Taylor instabilities, i.e., bubble/spike instabilities that drive small-scale mixing and reaction, and (3) Kelvin/Helmholtz instabilities that arise from radial entrainment of cool air due to the buoyant vertical acceleration of the pool fire.^{5}

In the LES regime, the production of kinetic energy and subgrid modeling approaches for buoyancy instabilities are extremely limited.^{42} Therefore, in most cases, these fundamental instabilities must be resolved in order to capture high-level pool fire phenomena such as fingering, puffing, and general core collapse.^{5,31} For example in a hybrid-RANS-LES simulation study of a non-reacting helium plume,^{10} mesh resolution for a low-dissipation second-order spatially accurate scheme required O(1) cm resolution to adequately capture centerline velocity core collapse. Therefore, as pool fire diameters increase to O(5) m and beyond, in the absence of modeling advances, modern large-scale compute resources must be exercised.

The requisite modeling of participating media radiation coupling physics for fire applications supports the exploration of many techniques, e.g., discrete transfer, discrete ordinates, spherical harmonics, Monte Carlo-based approaches, etc. (see Tencer and Howell^{43} for an overview of common techniques used for participating media radiation couplings or Jensen *et al.*^{44} for a fire use case comparison). In this study, and as described in Sec. II C, the discrete ordinate method^{45} represents the approach that we have deployed to solve the radiative transport equation (RTE). This approach allows for a good compromise between accuracy and, given our generalized unstructured implicit solver paradigm, code implementation simplicity. The discretization approach commonly used for complex geometries is the finite volume method^{46,47} with diamond mean flux,^{48} exponential,^{49} or upwinded^{50} advection operators (see the study by Joseph *et al.*^{51} for a comparison of various techniques). Upwind-based, finite-volume methods can also include gradient reconstruction and limiting techniques^{52} to remove unwanted oscillatory solution realizations. However, finite element-based approaches based on the Galerkin method are also noted, especially in complex geometries.^{53,54} Residual stabilization methods in the finite element context have also been demonstrated as a viable approach^{55,56} and provide a consistent approach whereby a fine-scale residual is evaluated to support the RTE stabilization strategy. In this study, we develop a consistent residual-based, finite-volume technique using a control-volume, finite-element method (CVFEM) to avoid *ad hoc*-based advection operators and remove the usage of limiter-based upwind methods to control monotonicity of integrated quantities such as scalar flux.

The Sierra/Fuego simulation tool^{6} that is exercised in this study was developed to support the capability-class computing, fire-physics use case. Fuego is part of a suite of numerical simulation tools used to address the abnormal thermal environment that is characterized as an accident scenario requiring thermal response modeling for objects exposed to fire. The Fuego code base is built under the Sierra Toolkit^{57} and Trilinos^{58} infrastructure that were first established in the extreme-scale computing regime using our open-source sister code, Nalu.^{29,59} The generalized unstructured, design-order, low-dissipation control-volume, finite-element method (CVFEM) outlined in the study by Domino *et al.*^{25,28} is implemented within Fuego. The Fuego simulation tool also holds a rich pedigree of code verification^{60} and is the tool that has been deployed to previous pool fire simulation studies,^{12} liquid natural gas^{61} accident scenarios, and other high-consequence simulation studies including nuclear reactors^{62} and pathogen transport using Eulerian/point-Lagrangian multi-physics coupling methods.^{63,64} As noted, the combustion modeling approach ranges from the already cited EDC usage to our flamelet-based approach outlined in this paper.

In this paper, our primary objective is to exercise a state-of-the-art combustion modeling methodology in the context of a large-eddy simulation turbulence closure paradigm to establish the credibility of this model suite's ability to capture key fire dynamics (puffing and fingering) and to assess the quality of our primary quantity of interest, radiative heat flux. We view the inclusion of a solution sensitivity methodology, i.e., the quantification of numerical sensitivity to mesh and time step refinement, to be a critical step in providing mesh resolution requirements for this class of buoyantly driven flows.

We first outline the numerical and math model approach for high-fidelity pool fire modeling, see Sec. II. This primary section includes the detailed variable-density, low-Mach number turbulent reacting approach and captures our discretization and coupling methodology. Where code verification exposition exists for the low-Mach scheme, references are provided. For the novel participating media radiation solver presented, detailed code verification is presented to establish the design-order accuracy of our proposed method. The full combustion model, which is based on a generalized heat loss flamelet approach, is outlined in Sec. III. Soot and radiative source terms are described in this section, in addition to tabulated approaches for these high-sooting hydrocarbon fires. The validation study (with mesh refinement sensitivity) of a 5 m pool fire, which is modeled as a heptane fuel, is provided in Sec. IV. Finally, Sec. V captures overall conclusions and paths to explore related to exascale computing initiatives and turbulence model improvement.

## II. NUMERICAL AND MATH MODEL OVERVIEW

As noted above, the fire simulation tool developed and deployed in this study is the National Nuclear Security Administration, Advanced Simulation and Computing massively parallel, generalized unstructured, turbulent low-Mach code base.^{6} This tool, which is named SIERRA/Fuego, is designed for low-Mach object-in-fire scenarios that include fuel-, propellant-, and composite-fires with support for fire suppression use cases. The reacting flow, variable-density, low-Mach number equations are coupled to a participating media radiation capability. Note that a low-Mach physics set is simply characterized as flow whereby the acoustic pressure waves have been filtered^{65} as they are generally not of interest to the thermal exposure question posed. Therefore, the equation of state for density evaluation is a function of the constant-in-space, variable-in-time (for a closed system) thermodynamic pressure. Subgrid stress model closure is obtained by a suite of large-eddy simulation models including the Smagorinsky model (with and without a dynamic procedure^{66}) and, in this simulation study, the one-equation *k _{sgs}* static coefficient approach.

^{67}For more information on this tool, the reader is referred to previous works

^{6,68}and the extensive approximate projection numerical methodology developed in both Fuego and its open-source sibling code, Nalu.

^{25,28,59}

### A. Low-Mach Reacting Flow Equation

The details of the mathematical model suite begin with the integral form of the Favre-filtered, variable-density continuity and momentum equations,

and

where the subgrid scale turbulent stress $\tau ijsgs$ is defined as

The stress is defined by,

and the traceless rate-of-strain tensor defined as follows:

For reacting flow, the integral forms of the Favre-filtered static enthalpy energy equation $h\u0303$, mixture fraction $Z\u0303$, source enthalpy $h\u0303c$, soot moles per mass $Ns\u0303$, and soot mass fraction $Ys\u0303$ used are as follows:

where $Sc$ and $Pr$ are the Schmidt and Prandtl numbers, respectively. For each equation, the subgrid scale turbulent diffusive flux vector $\tau \varphi ,jsgs$ is defined, for a representative variable $\varphi $, as

A gradient-diffusion closure model is used for $\tau \varphi ,jsgs$ in the context of isotropic, eddy-viscosity large-eddy simulation closure. A detailed description of the source terms is provided in Sec. III B. Note that the mixture fraction equation has an additional source term, $S\u0307Z$, due to movement of carbon from the gas phase to the soot mass.^{23,69,70} The mixture fraction is therefore a “gas phase” mixture fraction, tracking fuel from the inlet that remains gaseous while the soot mass fraction tracks the remainder. The term limits the growth of soot in the high-residence-time regions of the fire where otherwise the soot mass can grow unphysically large.

We introduce the source enthalpy, *h _{c}*, for use with the flamelet library approach to provide a reference enthalpy state based on the source oxidizer and fuel enthalpy streams. Radiative and wall heat losses are significant in fires, and this source enthalpy provides a link to the source stream state that is the flamelet reference state described in Sec. III C. An important source of variation in the fuel source enthalpy, for example, is the variation in the pool temperature and the amount of heat absorbed by the pool when the fuel vapor pressure is below atmospheric pressure. The source enthalpy is a conserved scalar that acts to increase the dimensionality of the flamelet concept to address the need for varying fuel and oxidizer source enthalpies.

### B. Discretization and coupling

For spatial discretization, the code base exercises the generalized, unstructured, control-volume, finite-element method, which has recently been verified and authenticated for large-eddy simulation usage.^{25} In this method, which combines the strengths of finite volume and finite element approaches, a dual mesh is constructed within the element.^{71} The formulation supports both low-order dual volume definitions for the hexahedral, tetrahedral, wedge, and pyramid element types and allows for higher-order elements such as a quadratic basis for the hexahedral 27-node element.^{72} Dual volume surface quadrature points, which evaluate integrated-by-parts terms such as advection and diffusion contributions, are defined on the sub-control surfaces while volumetric quadrature points, which evaluate volumetric terms such as time and source contributions, are defined at each sub-control volume. For more details on the method exercised in this study, the reader is referred to the recent works of Domino *et al.*^{25,28}

For the time-split, time-accurate, large-eddy simulation, regime-coupling approach exercised in this paper, we use an equal-order, pressure-stabilized, approximate pressure projection-based scheme.^{68} Variable density code verification of the methodology can be found in the study by Domino *et al.*^{25,60}

### C. Thermal radiation

For applications that include combustion and thermal radiation due to soot- and gas-phase energy participation, the divergence of the radiative flux is required for the gas-phase energy equation. The spatial variation of the radiative intensity corresponding to a given direction and at a given wavelength within a radiatively participating material, *I*(*s*), is governed by the Boltzmann transport equation, here shown with isotropic scattering activated,

Above, *α _{a}* is the absorption coefficient,

*α*is the scattering coefficient,

_{s}*I*(

*s*) is the intensity along the direction

*s*,

_{i}*T*is the temperature, and the scalar flux is

*G*. The black body radiation,

*I*, is defined by $\sigma T4\pi $. Note that for situations in which the scattering coefficient is zero, which is a good approximation for fire applications, and for configurations that include nonreflecting boundaries, the radiative transport equation reduces to a set of linear, decoupled equations for each intensity.

_{b}As noted in the introduction, the method of discrete ordinates^{45} is used to solve the radiative transport equation system. The symmetric quadrature set of Thurgood *et al.*^{73} is used for the quadrature set order, *N*. This quadrature set choice results in $8N2$ ordinate directions solved during each RTE iteration.

The novel discretization approach used in this study can be derived under the variational multi-scale (VMS) framework,^{74} however, to be cast within the control-volume, finite-element methodology. In the VMS framework, the degree of freedom and test function are decomposed in terms of its resolved, also known as, coarse-, and fine-scale feature, e.g., $I+I\u2032$. This decomposition is substituted into the RTE equation where the resolved and fine-scale terms are collected. Grouping resolved and fine-scale terms, below shown for the coarse-scale statement, results in an equation that takes the form of a standard Galerkin contribution in addition to a fine structure for the intensity degree-of-freedom contribution,

where *w* is the test function, and integrals are over the full domain, Ω.

The first term in the second line of Eq. (9) is integrated by parts to yield the following form:

The integral over Γ represents exterior surfaces, or boundaries of the domain. The following ansatz, obtained from the fine-scale statement, introduces the classic stabilization parameter, *τ ^{h}*, and provides closure for the above fine-scale degree of freedom,

Here, *R*(*s*) is the fine-scale residual for the governing RTE without which only the un-stabilized Galerkin form resides.

A compact form of the equation is provided by defining a modified test function, $w\u0303$,

where $w\u0303$ is simply equal to,

When $\alpha =\u22121$, we have the above VMS derivation; for *α* = 1, Galerkin Least Squares is realized; finally for *α* = 0, which is used in this study, we have replicated the streamwise upwind Petrov–Galerkin (SUPG) methodology generally used for advection-dominated partial-differential systems. For any formulation other than VMS, the residual contribution at the boundaries of the domain is dropped (*β* = 0).

The full residual-based equation is placed in divergence form,

and split into its Galerkin and stabilized contributions,

Note that the first term in the above equation is integrated by parts,

Again, the usage of Γ provides emphasis that the contribution is a boundary condition, i.e., an exposed face. Therefore, the full VMS-based stabilized RTE equation is as follows:

#### 1. Definition of the CVFEM test function

Following the work of Martinez,^{75} the test function is chosen to be a piecewise-constant value within the control volume, *w* = *w _{I}*, and zero outside of this control volume (Heaviside). A key property of this function is that its gradient is a distribution of delta functions on the dual volume boundary,

where Γ_{I} is the boundary of control volume *I* and **n**_{I} is the outward normal on that boundary. Substituting this relationship into the residual equation provides the final form of the vertex-centered, finite-volume, stabilized RTE equation,

Given this equation, with the test function defined accordingly, either an edge-based or element-based scheme can be used as outlined in the non-conformal hybrid methodology to convert a Discontinuous Galerkin method to a vertex-centered suite of methods; see Domino.^{28} For *α* = 0 and *β* = 0, which are used in this study, a Streamwise Upwind Control Volume (SUCV) is obtained.^{76} A residual-based stabilization methodology, then within a finite element context, was first introduced in Burns^{56} to solve the RTE in this manner.

For the boundary condition implementation, when $sini$ is greater than zero (assuming the outwardly normal boundary standard), the interpolated intensity based on the surface nodal values is used. However, when $sini$ is less than zero, the intensity boundary condition based on a diffuse surface is used, i.e.,

Above, *τ* and *ε* are the transmissivity and emissivity, respectively, while the background reference temperature and surface temperature are $T\u221e4$ and $Tw4$, respectively. Finally, $Hr$ is the irradiation, or sum of incoming intensities to the surface. Since the original RTE equation was integrated by parts, a natural surface flux contribution is applied. For general reflecting surfaces, the set of all intensity equations are linked to the calculation of the surface irradiation.

#### 2. The form of *τ*^{h}

^{h}

The value of the stabilization parameter *τ* can take on a variety of forms. In general, the limiting forms of *τ* represented in the advection and diffusive limit are $\tau advh=h2$ and $\tau diffh=1(\alpha a+\alpha s)$, respectively. In this study, we follow the finite element method SUPG-based approach of Avila *et al.* that also provides good accuracy and stability in a simplified and easy to implement form,

In the above equation, the element length scale, *h* is determined by the dot product of the ordinate direction *s _{i}* and the area vector unit normal at the control volume surface integration point.

#### 3. Verification

A detailed code verification study is presented to establish the order-of-accuracy for the SUCV edge- and element-based advection operator using the residual-based stabilized approach outlined in this paper. The multi-physics verification test case represents a previously established thermal/PMR coupling in the thin limit, see the study by Domino *et al.*,^{60} that used our formerly deployed the SUPG methodology^{56} in the absence of radiative source terms. This verification test defines four nested concentric spheres of different radii that form two separate heat conduction regions separated by a single vacuum region, i.e., nonparticipating media radiation, of gap 0.25 m [see Fig. 1(a)]. The total radius of the outer region is 0.1 m. The internal thermal boundary condition links the computed irradiation from the PMR solve to the thermal response,

where *λ* is the thermal conductivity, *n _{j}* is the outward pointing surface normal, and

*G*is again the surface irradiation. In this three-dimensional configuration, the inner-most temperature is set to 300 K while the outer surface is 600 K. Interior emissivities are 0.8 and 0.6 for the inner and outer heat conduction surface, respectively. As noted, this test provides verification coverage for the advection operator and exercises the mesh transfer infrastructure between disparate meshes for each physics as provided by the Sierra Toolkit Transfer module.

In Fig. 1(b), the edge- and element-based vertex-centered approach is used for linear hexahedral and tetrahedral topologies (Hex8 and Tet4 element topologies, where 8 and 4 represent the number of nodes in the element) in addition to a 27-node quadratic hexahedral element type (Hex27). The *L*_{2} norm is based on the analytical and simulated temperature field within the outer thermal domain. For the low-order element topologies (Hex8 and Tet4), a Thurgood quadrature set using order 10 (*T*10) was used while for the quadratic CVFEM basis, a quadrature order of 22 (*T*22) was activated. Sub-optimal convergence due to quadrature error pollution was found when using the *T*10 quadrature order for the quadratic test case as illustrated in the figure. The usage of the lower quadrature order and its effect on convergence of the third-order spatially accurate quadratic basis CVFEM simulation [see Fig. 1(b)] demonstrates the dominance of quadrature error on the *L*_{2} norm. Therefore, design-order, i.e., a second-order slope for a linear basis and third-order slope for a quadratic basis, can only be established for this challenging, optically thin, multi-physics verification case using the proposed SUCV-based stabilization when quadrature errors are reduced below numerical discretization errors. Also noted in this simulation study, however not presented, is the requirement to refine the mesh while adapting to the spherical geometry fidelity to avoid errors that are dominated by faceting of the curved surface.

## III. COMBUSTION MODEL SUITE

To predict the evolution of the thermochemical state, a generalized flamelet based approach is employed. While the present approach has many similarities to traditional flamelet modeling, the importance of soot and radiation leads to some significant differences as detailed in Secs. III A–III C.

### A. Flamelets with heat losses

The majority of fuel oxidation rates in fires are substantially faster than the flow mixing rates, so that combustion is firmly in the nonpremixed combustion regime. The mixing process is then described using the mixture fraction, *Z*, that represents the fraction of the mixture originating in the fuel stream. The flow mixing rate is characterized using the scalar dissipation rate, $\chi =2DZ|\u2207Z|2$ where *D _{Z}* is the relevant diffusion coefficient for the mixture fraction (obtained from the viscosity with a Schmidt number of 0.9). Quantities needed by the simulation are referenced to the mixture fraction, the scalar dissipation rate, and the enthalpy because of the significance of radiation heat transfer in fires. Needed quantities include source terms for soot formation and oxidation, radiative source terms given in Sec. III B, and the density and viscosity that are involved in all of the transport equations in Sec. II A. This abstraction conveniently separates the fluid solution to the evolution of those variables that are relevant at fire evolution time scales, which are the fuel–air mixing time scales, from the faster fuel oxidation processes that are precomputed as described here.

To develop the above needed quantities, as well as any other description of the thermochemical state, flamelet equations are solved for the fuel–air mixing process using a chemical kinetic mechanism for heptane as a fuel surrogate.^{77} A series of adiabatic steady flamelets is obtained spanning the range of scalar dissipation rates from very low values important in fires (as low as $\chi st=10\u22123\u2009s\u22121$) up to $\chi st=50\u2009s\u22121$ near the extinction condition. Each of these flamelets is then evolved using the unsteady flamelet equations that include a heat loss term for the energy equation,

wherein *Y _{i}*,

*h*, $cp,i$, and

_{i}*ω*are the mass fraction, enthalpy, isobaric specific heat capacity, and net mass production rate of species

_{i}*i*, respectively,

*T*is the temperature,

*t*is the flamelet evolution time,

*ρ*is the mass density,

*c*is the isobaric specific heat capacity of the mixture, $T\u221e$ is a chosen cooled reference temperature, and

_{p}*H*is an arbitrary heat transfer coefficient, discussed further below.

Although the early and traditional unsteady flamelet formulations introduce a flamelet timescale that is determined based on the distance from the inlet/nozzle,^{78} termed a Lagrangian unsteady flamelet formulation,^{79} the Lagrangian approach is less appropriate in more elliptic buoyant flows. Our approach follows a method closer to Ihme and Pitch^{80} where unsteady flamelet solutions are carried out separately and provide appropriate source terms needed for full elliptic solves of the soot and enthalpy. Specifically, for each dissipation rate, the flamelet is evolved from the adiabatic condition until the stoichiometric temperature is within a few percent of the frozen, i.e., cold-mixing, stoichiometric temperature. The resulting solution set spans the full range of dissipation rates and possible enthalpies, from adiabatic to fully cooled. A range of example profiles for the enthalpy and the temperature is provided from the resulting flamelet solutions in Fig. 2 for the scalar dissipation rate $\chi st=0.1032\u2009s\u22121$. This set of solutions is then fit to a multidimensional table that maps filtered source terms to filtered evolved quantities by a process described in Sec. III C.

The heat loss in the flamelet equation need not be directly linked to a specific heat-loss mechanism (like radiative emissions), thereby providing a flexible approach. For this application, a heat loss that is linear in temperature is used since this approach provides a good range of values for the given fire conditions. For many hydrocarbon fuels, fuel-pyrolysis chemistry creates positive curvature in the temperature vs mixture fraction curve on the rich side of the flame as seen in Fig. 2; this curvature is reflected in the enthalpy profiles as well. If the heat loss term is made proportional to $T4$ as expected for radiative losses, the curvature in the enthalpy profiles would be more severe. As discussed in Sec. III C, the shape of the enthalpy profiles determines the effective correlation between a local flame state and the important stoichiometric flame state.

In the heat loss term, there is a product of the scalar dissipation rate with a coefficient, *H*. This product appears because the heat losses are designed to overcome the heat release to allow complete cooling at each dissipation rate. With mixing-limited heat-release rates, the flamelet heat-release rate is proportional to *χ* so that a cooling rate for a low dissipation rate will be insufficient for a high dissipation rate flame, hence the proportionality to *χ*. Conversely, the cooling required for a higher dissipation rate flamelet would completely overwhelm diffusive mixing for a flamelet with a small scalar dissipation rate leading to profiles of enthalpy heat-loss that would be sharply peaked around the stoichiometric point. The heat loss required to cool a flame is then somewhat greater than the flame heat release rate itself. The latter can be estimated with the dominant diffusive term of Eq. (25), $(\chi /2)\u22022T/\u2202Z2$, using an estimate that follows.^{81}

The change in the temperature gradient across the flame is proportional to $(Tmax\u2212T\u221e)Zst\u22121(1\u2212Zst)\u22121$, divided by a reaction-zone thickness, so an estimate of the heat release rate is proportional to these parameters. The parameter *H* is related to the flame temperature rise $(Tmax\u2212T\u221e)$, or equivalently the flame heat release, divided by a reaction zone thickness giving,

Because the heat loss must exceed the heat release rate by some margin, there is some empiricism in the particular value of *H* identified, but we have found that the value does not vary substantially over a range of hydrocarbon flames of interest to our group.

In Fig. 3, we visualize characteristics of cooling flamelets in a normalized flamelet evolution time through the maximum values of temperature and OH mass fraction, for two values of the dissipation rate and three values of *H*. This shows that flamelets are cooled until peak flamelet temperatures drop below 500 K, where cooling causes a collapse of the radical pool. This collapse of the radical pool is reflected in the maximum OH mass fraction, important for soot oxidation, approaching zero. Thus, for the greatest heat losses, the flamelet has undergone a cooling driven extinction where soot oxidation will be minimal, allowing smoke emissions. Figure 3 also shows $\gamma st$ over the normalized flamelet evolution time to demonstrate how choice of *H* determines the range of heat losses included in the flamelet table. Flamelet tables used in this study employed $H=107$ to encompass a large range of heat loss spanning more than $\gamma st\u2208[0,\u22122\u2009MJ/kg]$.

Beyond radiative extinction times, reactants diffusing into the flame are not consumed, and the stoichiometric value for the chemical enthalpy increases in time; this is associated with the diffusion of products away from the stoichiometric region and the diffusion of reactants toward the stoichiometric region. This trend follows an exponential approach to the frozen mixing enthalpy that is the same as the adiabatic flame enthalpy. Simultaneously during the cooling process, the temperature continues to decay, leading to a continued reduction in the sensible enthalpy. When we use the above source term with $Tmax(t)\u2212T\u221e$ varying with the current value, then the temperature reduction is linear in time as compared with an exponential decay that occurs if a fixed value of $Tmax(t=0)\u2212T\u221e$ is used in the denominator of the heat loss term of Eq. (24). This linear-in-time sensible enthalpy behavior allows the unsteady flamelets to expand to a broader range of enthalpy space where the late stage cooling is dominated by cooled reaction products rather than reactants.

As the cooling rate is reduced with smaller *H* coefficients, more reactant mixing occurs as shown in the *γ _{st}* curves for smaller

*H*in Fig. 3. The lowest value of

*H*leads to 25 to 35% reactants in the cooled products while the largest value of

*H*is dominated by cooled products at the greatest extent of cooling. We also observe that

*γ*is beginning to evolve toward a positive slope for the smallest

_{st}*H*, which would lead to a non-monotonic enthalpy vs cooling behavior. We note that we do not have any direct evidence that one value of

*H*is better than another, but the LES mixing itself directly predicts the mixing of reactants. The role of the flamelet models is to generate a range of states that might be needed in general scenarios, for example where turbulent mixing rates are low and heat losses are large. With some flexibility in selecting

*H*, this shows that it is possible to vary the extent to which the reactant/product ratio varies with the greatest degree of cooling.

### B. Soot and radiative source terms

The flamelet solutions developed in Sec. III A are used to generate values needed for the LES, including parameters like the density and source terms for evolving variables. Source terms are needed for the enthalpy and for soot, both of which evolve over slower time scales comparable to resolved fluid mixing.

Soot is modeled using a two equation model based on the work of Aksit and Moss.^{82} Transport equations modeling the soot moles per mass and the soot mass fraction in the flame are given in Eq. (6), where $Sc=Pr=0.9$ and $Scs=200$. Diffusion of soot by Brownian motion implies a very large Schmidt number for soot, e.g., $Sc=290$ for 10 nm spherical particles;^{83} the value of 200 is chosen to be sufficiently large where only turbulent diffusion is important. One modification is needed for the source terms in large turbulent fires: in the core regions of the fire, rich, soot-laden mixtures can absorb radiant energy and remain quite hot. Because soot formation in these two equation models follows a simple Arrhenius rate, elevated temperatures lead to unphysically large soot production rates. It is known that soot growth follows the so-called bell curve (e.g., Frenklach *et al.*^{84} and Alexiou and Williams^{85}) where molecular decomposition inhibits soot growth at high temperatures. To mimic this behavior, we introduce one additional term, $\omega \u0307R$, that smoothly adjusts the soot nucleation and surface growth to a negative value above temperatures of *T _{c}*; the rate of this transition is given by an effective decomposition activation temperature,

*T*.

_{R}In fires, an appreciable fraction of the fuel carbon can be associated with the soot. This has a significant impact on the fuel-consumption process. To address this, soot mass formation is linked to the mixture fraction by providing the same source terms of opposite sign. Without this change and the introduction of $\omega \u0307R$, unphysical levels of soot are predicted in large-scale fires.

The source terms for soot moles per mass, soot mass fraction, the mixture fraction, and the enthalpy from Eq. (6) are

where $\alpha gas$ is determined from the mixture-averaged absorption coefficients of $CO,\u2009CO2,\u2009CH4$, and $H2O$, curve-fit from RADCAL data.^{86,87} Closures for soot production terms, e.g., $S\u0307coagulationYs\u03031/6N\u0303s11/6$ are neglected in this study; however, the subgrid fluctuations have been found to be important in smaller scale pool fires.^{37}

The mass-fraction-normalized absorption coefficient of soot is taken to be a linear function of temperature, $\alpha soot=max(\alpha 0+k\alpha T,0)\rho /\rho s$ with $\alpha 0=\u2212375\u2009000\u2009m\u22121,\u2009k\alpha =1735\u2009m\u22121\u2009K\u22121$, and $\rho s=1800$ kg m^{−3}.^{88} The radiation sources are $qsr=\alpha sootT4$ and $qgasr=\alpha gasT4$. The soot-related rates follow^{82} with the addition of the new $\omega \u0307R$ term

where $TR=16\u2009000\u2009K,\u2009Tc=2050\u2009K,\u2009AN=54\u2009s\u22121,\u2009TN=21\u2009100\u2009K,ASG=12\u2009000\u2009s\u22121,\u2009TSG=12\u2009000\u2009K,\u2009AO2=500\u2009s\u22121,\u2009TO2=20\u2009000\u2009K,AOH=4.2325\u2009s\u22121$, and $\rho s=1800\u2009kg\u2009m\u22123$. Above, $vd=61/3(\pi \rho sNA)\u22121/3$ converts volume to an effective particle diameter together with the $Ys\u03032/3N\u0303s1/3$ to give a surface area dependence for the soot surface growth and oxidation terms; *N _{A}* is Avogadro's number. The $O2$ oxidation term is neglected by Aksit and Moss,

^{82}and model parameters from Leung

*et al.*

^{89}are used to include its effect. The nucleation and coagulation source terms are in terms of moles per mass of mixture, and

*W*= 144 kg kmol

_{s}^{−1}is the assumed nucleate mass to describe the nucleation contribution to the soot mass fraction. The overline in each of the above expressions (including the radiative source terms) denotes an averaging of quantities over subgrid fluctuations using a presumed probability-density function (PDF) approach that is precomputed as described in Sec. III C.

The parameters given above are generally associated with ethylene flames, and JP-8 is expected to produce soot significantly faster than ethylene. However with the slow fuel–air mixing at these scales (see Sec. IV), the soot levels, optical thicknesses, and radiative heat fluxes are found to be in line with measurements. We expect that fires with length scales comparable to tens of centimeters will have a much greater sensitivity to the soot production rates. Soot mass fractions in a large-scale pool fire can range up to several percent on average, with much higher instantaneous values (see Fig. 11). Due to the large absorption coefficient of soot, this fire is firmly in the optically thick regime where participating media radiation treatments introduced here are important for predicting the environment.

### C. Tabulation approach

The LES filtering process leads to unresolved mixing that is important in high Reynolds number flows. This mixing is accounted for in all flamelet-derived quantities using a presumed-pdf approach (e.g., Peters^{90}) including $\rho \xaf$ and all of the source terms in Sec. III B; the overline extent shows the terms accounted for in the filtering process. Terms with an overline are evaluated using

where $PZ\chi \gamma $ is the joint PDF of the sub-filter fluctuations of $\varphi $ in the $Z\u2212\chi \u2212\gamma $ space. Presuming statistical independence of the *Z*, *χ*, and *γ* fluctuations, and letting the dissipation rate and enthalpy defect PDFs be modeled as delta functions, yields the more tractable expression,

in which *P _{Z}* is the mixture fraction PDF. The convolution integral above is preprocessed and tabulated for the desired range of mixture fraction and scalar variance values. The Favre-filtered property $\varphi \u0303$ could then be evaluated directly from the Favre-averaged inputs as in the left-hand side of Eq. (38) using a multidimensional, unstructured interpolation or regression technique. In order to avoid the complexity of high-dimensional, unstructured interpolation and the inaccuracy of regression approaches, we exploit the nature of the flamelet-library generation process to enable the calculation of $\varphi \u0303$ with efficient, accurate, structured interpolation techniques.

Results of this filtering operation are organized into a lookup table for efficient evaluation of tabulated properties such as density and soot source terms. For laminar Burke–Schumann or equilibrium chemistry models, these tables are one-dimensional in nature, and piecewise polynomial interpolation is straightforward for tabulation. Here, the flamelet library exists in four dimensions, $Z\u0303,\u2009Z\u20332\u0303,\u2009\chi \u0303$, and $\gamma \u0303$, that do not immediately fill a structured mesh as suggested by the shape of the enthalpy-defect profiles in Fig. 2. While there may be promise in unstructured regression techniques such as multivariate adaptive regression splines and artificial neural networks, we utilize a specialized technique that retains accuracy and leverages structured interpolation methods.

To enable the use of structured lookup tables, the $\chi \u0303,\u2009\gamma \u0303$, and $Z\u20332\u0303$ input values are converted to reference or scaled values $\chi st,\u2009\gamma st,\u2009Z\u20332s$, over which structured interpolants of properties can be built. These use the correlation inherent in the flamelet solution for the enthalpy defect and the scalar dissipation rate. For fires, the scalar dissipation rates are small (very far from extinction), and the results are weakly dependent on the profile of the scalar dissipation rate, but the correlation between the enthalpy defect and the mixture fraction depends strongly on the heat loss term introduced in Sec. III A. Generally, let $\eta \u0303$ represent the vector of filtered inputs acquired directly from evolved quantities in the simulation, and let $\eta o$ be a vector of input variables over which the flamelet data form a structured table. Conversion of filtered inputs to structured inputs can be accomplished by tabulating the filtered inputs, $\eta \u0303table(\eta o)$, and solving a minimization problem to identify the structured inputs that most closely reconstruct the filtered values:

In general, treating table lookup as a minimization problem atop a structured interpolation table is likely not feasible if multiple dimensions must be ‘inverted’ in this manner. Due to the use of a presumed form of the dissipation rate in building the laminar flamelet table, $\chi (Z)$, the inversion of $\chi \u0303$ to compute the corresponding stoichiometric value, $\chi st$, can be done explicitly. Similarly, an explicit conversion exists for the scaled scalar variance, $Zs\u20332=Z\u20332\u0303/Z\u0303(1\u2212Z\u0303)$.

Simultaneously including both finite strain and heat loss effects, however, leads to flamelet libraries where *γ* depends nonlinearly on *Z* and *χ*, precluding explicit input conversion and requiring the minimization approach to ensure the consistency of observed and tabulated enthalpy defects. We refer to this approach as consistent enthalpy reconstruction (CER), a technique needed for generalized heat loss models with finite strain. The minimization problem is now one-dimensional in the $\gamma st$ coordinate,

where the $Zs\u20332$ and $\chi st$ values have been converted explicitly from the filtered input variables. We solve this equation for consistent values of $\gamma st$ with a secant method modified to account for finite table boundaries, and we fall back to Brent's method if the secant solver fails to converge quickly. We observe that typically between three to six additional evaluations of the structured interpolation table are needed to identify the consistent value of $\gamma st$, which can be computed just once for a given set of filtered inputs and be reused over many properties, minimizing the additional cost.

The choice of the stoichiometric reference state is arbitrary, although it does facilitate matters when presumed forms of the dissipation rate or enthalpy defect are employed. The solution of the minimization problem (40) is not a requirement intrinsic to the heat loss approach. It is entirely a consequence of using structured, tensor product interpolation routines to evaluate properties that are defined on an unstructured dataset. At this point, the feasibility of unstructured, high-dimensional interpolation methods with at least a cubic polynomial basis or regression methods remains unclear.

## IV. 5 m POOL FIRE VALIDATION

In general, jet fuel is a complex mixture of underlying components including normal and branched alkanes, cyclo-alkanes, aromatics such as octane and dodecane, and bi-cyclic organics such as decalin.^{91} In this case study, a 5-m-diameter, *D*, heptane pool fire is used to model JP-8 fuel. For this configuration, the radiative heat flux at the pool surface has been measured and first reported as a RANS-based validation study by Domino *et al.*^{92} In this section, we define the computational domain and discuss the fire dynamics predicted as a function of computational mesh spacing. The study showcases the LES-based, unsteady flamelet approach with the proposed soot model in the context of our current gray participating media assumption. We also discuss measurements of a slightly smaller 2 m diameter JP-8 pool fire measured by Murphy and Shaddix^{93} due to the similarities between the mid- and large-scale pool fires.

### A. Domain definition and modeling approach

The experimental configuration is modeled as a cylindrical domain 8*D* in diameter and 10*D* in height. For the coarsest mesh, the inner-pool surface spacing was 0.1 m and transitioned to 0.05 m at the pool lip. An inner ring of refinement (from 1*D* to 1.5*D*) allowed for transitioning from the pool lip, i.e., the ground-level liquid/air interface, to the far-field entrainment boundary. The coarse mesh transitioned radially to 0.5 m mesh spacing at the entrainment boundary. The top domain mesh spacing at 10*D* above the pool was 0.5 m while the first node off of the pool surface was 0.1 m. The unstructured baseline mesh consisted of approximately two million linear hexahedral elements. The two subsequent meshes included in this refinement study were uniform mesh spacing reductions (each mesh resulting in a factor of two smaller mesh spacing) that provided element counts of approximately 16 and 128 million linear hexahedral elements, respectively. Therefore, the 128 million element mesh provides a pool lip mesh spacing of 1.25 cm. Due to the re-meshing of the domain, pool curvature was captured during the mesh refinement process. The finest resolution simulations were run on the Eclipse Sandia National Laboratory institutional cluster using 128 nodes (2.1 GHz processors; nodes with dual sockets, 18 cores each; Intel Omni-Path high-speed interconnect) for a total of 4608 Message Passing Interface ranks. Approximately ten 48-h sessions were run for the time-mean plots provided in this section. Figure 4 provides a wireframe overview of the most-coarse mesh used in this study (5.0 cm nominal resolution at the pool lip).

For turbulence modeling, the one-equation *k _{sgs}* model was used.

^{67}Low-dissipation CVFEM-based advection was activated following the work of Domino

*et al.*

^{25}while using a second-order time accurate BDF2 time integrator with a time step of near unity Courant number. For scalar transport, a Van-Leer MUSCL-based stabilization approach (with limiter) was used.

^{6}The radiative transport equation was solved using edge-based SUCV stabilization with a quadrature order of four (

*T*4 = 128 ordinate directions), which has been shown to be sufficient for near-fire radiative heat flux calculations.

^{92}For the fire application space where scattering effects are small, the RTE uses a scattering coefficient of zero.

The mass evaporation rate is generally a function of pool fire shape, size, and fuel type^{94} with devolatilized gas-phase fuel components varying over the time frame of the inventory (or total volume of liquid fuel) combustion event.^{91} In this study, a constant value of 0.053 $kg/s/m2$ was specified (temperature at 488 K) over the full simulation time^{95} while no pool liquid modeling was included as found in Fukumoto *et al.*^{96} At the pool surface, the transmissivity is specified to be unity allowing for a weak boundary condition that specifies the intensity at the pool surface to be based on the boiling pool temperature. Wall and open boundaries for the RTE also apply a weak boundary value corresponding to a far-field temperature of 25 °C using an emissivity of unity. Simulations were run for roughly 30 s of physical time with statistics captured over approximately ten puffing cycles.

### B. Fire dynamics

To showcase the effect of mesh resolution on the fingering pool fire dynamics, see Fig. 5. In this context, fingering instabilities, which are given the name due to the fingerlike structures seen at the base of the fire, are related to the interaction between the radial entrainment toward the center of the fire and Rayleigh/Taylor instabilities. In this figure, the volume-rendered temperature field is shown with specific emphasis on how the richness of the instability is captured with increasing mesh resolution.

In Fig. 6, an evolution of one puffing cycle is presented for the most refined mesh simulation. In this sequence, we can see the influence of both the Rayleigh/Taylor instability, i.e., bubble-spike behavior, the Kelvin/Helmholz behavior of fast entraining fluid toward the center of the domain, and the baroclinic torque term ($1\rho 2\u2207\rho \xd7\u2207p$), which is responsible for the large-scale vortical structure rotating the flame sheet. Slight spillage of the heavy fuel is noted past the inflow surface due to a perfectly flat interface between the pool and surrounding entrainment surface. The evolution of the primary vortex ring structure is noted along with stretching and tilting of this foundational intermittent structure as it accelerates vertically.

In Fig. 7, the temperature time series at a centerline point (two diameters in elevation) is provided in addition to the power spectral density (PSD) for the mixture fraction. The peak frequency noted in the PSD is approximately 0.6 Hz, which compares well to the correlation of Cetegen and Ahmed,^{97} $1.5D12$. We note that uncertainty in the correlation is not provided although the authors report that “the agreement of the data is truly remarkable.”

### C. Pool fire validation and mesh sensitivity

In Fig. 8, the time-mean, Favre-averaged, centerline temperature and vertical velocity component for the three mesh resolutions are provided. For the centerline temperature plot, the expected centerline maximum temperature height based on experimental correlations^{41} is also shown. The correlation-based approach, which requires specification, predicts a plume height of 3.44*D* as compared to the simulation, 3.08*D*. However, with refinement the trend is toward the correlation prediction. Note that the correlation does not provide precise temperature values, rather only a location at which the peak temperature is likely to occur. For the purposes of this study, uncertainties in flame height are obtained by simple variation in heat release rate. Specifically, for a ±10% variation in the heat release rate, flame heights range between 3.09 and 3.8*D*. Favre-averaged vertical velocity predictions at the centerline location are 25 m/s at the normalized height *h*/*D* of 3.85, which provides for the large Reynolds number (nearly 600 000 based on pool diameter) typically noted for large-scale pool fires. With refinement, both the temperature and velocity fields show evidence of slightly reduced dissipation/mixing. Less-resolved simulations have temperatures that peak earlier and velocities that peak lower; this is also reflected in the radial profiles below.

For quantitative comparisons to experimental data, Fig. 9 compares the radial profile of the radiative heat flux to the pool surface between simulation and experiment for the three mesh resolutions run in this study in addition to the baseline RANS results.^{92} The experiment measured a peak surface radiative flux of approximately 56 $kW/m2$ occurring at a *r*/*D* of 0.30. In this figure, experimental error bars are provided that represent a uniform 5% error in the reported measured value. This uncertainty is associated with a Sandia-developed heat flux gauge^{98} that has been well exercised in other experimental studies.^{99} Due to the large-scale size of the pool fire, the experiment was conducted outside. Uncertainties associated with quiescent to low-speed crossflow, which can substantially effect pool mass devolitilization,^{99} are factors that have not been quantified in this particular data set and again highlight the significant challenge in conducting large-scale outdoor fire experiments. In general, we note reduced error between the predicted and experimental peak value as the mesh resolution increases. Specifically, the most refined mesh predicts a peak surface radiative heat flux at an *r*/*D* of 0.32 with a nominal value of approximately 50 $kW/m2$. We believe that our current LES-based results are more accurate than our hybrid RANS study found in Domino *et al.*^{92} In our previous study using first-generation models that used extremely coarse meshes (nominal pool mesh spacing of approximately 13 cm), the maximum radiative heat flux magnitude was predicted to be approximately 78 $kW/m2$, a value closer to the most coarse mesh LES result predicting a value of approximately 68 $kW/m2$. Moreover, the RANS approach yielded a lifted flame, i.e., detached from the pool surface, that resulted in a reduction in radiative flux beyond the $rD$ value of 0.5. Finally, the plume dynamics captured in this baseline RANS activity, which is only noted for historical context, resulted in further outward shifting of the peak magnitude. Our current results, which now includes an order of magnitude lower mesh spacing, demonstrate that lower mesh resolutions predict a higher peak radiative heat flux, which from a safety perspective, provides a more conservative result as most production simulation studies are expected to be under-resolved. We hypothesize that the reduced radiative heat flux with higher-resolution meshes, which are approaching the experimental measurement, is due to capturing of the intermittent pool fire dynamics, such as puffing and local radiative cooling of large vortical structures.

#### 1. Radial profiles as a function of mesh resolution

In this subsection, we provide radial plots for temperature, soot, and radiative flux at various pool heights to capture the effect of mesh refinement on non-experimentally measured quantities of interest to possibly serve as future benchmarks for the 5 m pool fire. Specifically in Fig. 10, radial temperature values at varying heights for the three mesh resolutions are shown to illustrate mesh sensitivities. The plots provided again demonstrate that, with the lowest mesh resolution, the temperature profiles peak at smaller *h*/*D*, while there is reasonable agreement with the mesh refinements indicated by 2.5 and 1.25 cm base resolution. We suggest that lower resolutions lead to somewhat greater dissipative mixing, especially low in the fire where length scales are the smallest, and also give the reduced flame heights evident in Fig. 8.

For radial soot compositions at varying heights (as a function of mesh resolution), see Fig. 11. For the two more refined meshes, the predicted peak mass fractions correspond to soot volume fractions in the range of 2 to 4 ppm (greatest peak values near h/D = 2). These values are in reasonable agreement with measurements of Murphy and Shaddix^{93} in smaller, 2-m-diameter JP-8 pool fire where most measurements are clustered near 1 ppm with some toward 2 ppm. For the lowest resolution, soot predictions using the flamelet model are substantially higher in the lower regions of the flame, being roughly three times greater than the other two resolutions. The relatively strong dependence on the resolution here suggests that the source terms require somewhat greater resolution than is often employed. The relatively high soot mass fractions at the highest normalized vertical distance for the most finely resolved simulation occur because the longer flame length also allows for more rich pockets at that height near its core collapse.

Finally, in Fig. 12, we provide the magnitude of the net radiative flux vector at varying heights for the three mesh resolutions. The net radiative flux vector, $qjr$, in the discrete ordinates setting is computed using a set of ordinate direction weights, *w _{k}*, the ordinate directions, $sjk$, and the local intensity, $Ik$, and is given by the expression, $\u2211kwksjkIk$. As seen from the plot, the net vector magnitude is near zero at the centerline (except for a small vertical component) but provides our simulation-based representation of the radiative flux outside the flame region. In the lowest mesh resolution with the faster development of the temperature and the larger soot mass fraction, the radiative flux is somewhat greater, although the magnitude of the difference is not as large as might be expected. Results indicate that near to the ground level, where mesh refinement is greatest, the 2.5 and 1.25 cm results are very similar to those already shown in Fig. 9.

## V. CONCLUSIONS

In this paper, we have outlined a high-fidelity, large-eddy-simulation modeling and simulation paradigm to support the evaluation of the ability to capture salient pool fire dynamics noted in large-scale experiments. Moreover, our unsteady flamelet model with generalized heat loss was outlined and showcased for a highly sooting, large-scale, 5-m pool fire configuration under a gray assumption. To support multi-physics coupling, a novel, residual-based stabilization approach exercised within a control-volume, finite-element scheme and a two-state (left-right) edge-based, vertex-centered discretization for the RTE in the presence of non-scattering, participating-media radiation was presented and verified to be design-order for low-order hexahedral and tetrahedral element topologies in addition to the 27-node quadratic hexahedral element type.

With appropriate mesh resolution of nearly 1.25 cm, fire dynamics such as fingering and puffing are captured accurately. Specifically, the 5 m puffing frequency was well predicted as compared to the experimental data correlation-based approach of Cetegen and Ahmed^{97} (0.6 Hz vs nominal 0.67 Hz). Peak temperature centerline values were predicted at a plume height of 3.08*D* and compared well to the range of correlation-based heights. Moreover, and perhaps most importantly to our desired quantity of interest prediction, the model suite predicted realistic radiative heat flux values at a the pool surface as compared to experimentally measured values (again, approximately 10%). We note that the extensive usage of mesh and time step refinement techniques deployed in this paper is critical, without which no statement as to numerical sensitivity to end quantities of interest can be provided. However, due to the challenging requirement to capture key hydrodynamic instabilities such as Rayleigh/Taylor and Kelvin/Helmholz phenomena, continued capability-class computing investments must be made to allow for efficient solver strategies. Moreover, advanced turbulence model closure that relaxes the resolution requirement represents a rich investment application space.

As compared to previous RANS-based simulations using an EDC combustion model suite, increased predictivity was noted over the reduced fidelity turbulent and combustion modeling paradigm. Although solution sensitivity to underlying mesh resolution was exercised, mesh convergence was not perfectly established, even at the nominal 1.25 cm mesh resolution. However, key elements of fire dynamics such as core collapse were adequately predicted at all mesh resolutions. In fact, mesh refinement studies are planned using the capability-class computing resources available to the laboratory with expected mesh counts to surpass the one-billion element mark.

Finally, due to the conservative nature of predicting higher radiative heat fluxes at under-resolved simulations, production simulations on O(5000) compute cores represent a viable simulation size to carry out detailed fire analysis for large-scale fire use cases of interest to the laboratory. A fourth mesh resolution and relaxation of the gray assumption represent future activities.

## AUTHORS' CONTRIBUTIONS

S.P.D. is the lead author and Principal Investigator for the Verification and Validation portfolio that drove this research. J.H. is the PI for the fire Physics and Engineering Modeling portfolio that drove combustion and soot modeling improvements. R.K. and M.K. supported the code implementation for the soot and multidimensional table lookup code infrastructure. All authors contributed to the writing of this manuscript.

## ACKNOWLEDGMENTS

Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under Contract No. DE-NA0003525. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government, No. SAND2021-8617 J. This research was supported by the Advanced Simulation and Computing portfolio at Sandia National Laboratories. The authors thank graduate intern Elizabeth Armstrong and her advisor Professor James Sutherland (University of Utah) for collaborative work in support of the consistent enthalpy reconstruction formulation. Finally, this work is dedicated to the late Dr. Sheldon Tieszen whose vision and dream to advance Sandia's fire program has been our honor to serve.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article