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.

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, FDS19 or the Open-FOAM-based tool, Fire-FOAM20,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 studies25 that include structural uncertainty quantification26,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 velocimetry11 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 temperatures2,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 Howell43 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 method45 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 method46,47 with diamond mean flux,48 exponential,49 or upwinded50 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 techniques52 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 approach55,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 tool6 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 Toolkit57 and Trilinos58 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 verification60 and is the tool that has been deployed to previous pool fire simulation studies,12 liquid natural gas61 accident scenarios, and other high-consequence simulation studies including nuclear reactors62 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.

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 filtered65 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 procedure66) and, in this simulation study, the one-equation ksgs static coefficient approach.67 For more information on this tool, the reader is referred to previous works6,68 and the extensive approximate projection numerical methodology developed in both Fuego and its open-source sibling code, Nalu.25,28,59

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

ρ¯tdV+ρ¯ũjnjdS=0
(1)

and

ρ¯ũitdV+ρ¯ũiũjnjdS=σ̃ijnjdSτijsgsnjdS+(ρ¯ρ°)gidV,
(2)

where the subgrid scale turbulent stress τijsgs is defined as

τijsgsρ¯(uiuj̃ũiũj).
(3)

The stress is defined by,

σij=2μ¯S̃ij*P¯δij,
(4)

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

S̃ij*=S̃ij13δijS̃kk=S̃ij13ũkxkδij.
(5)

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

ρ¯Z̃tdV+ρ¯ũjZ̃njdS=(μ¯ScZ̃xjτZ,jsgs)njdS+ṠZdV,ρ¯h̃tdV+ρ¯ũjh̃njdS=(μ¯Prh̃xjτh,jsgs)njdSṠhdV,ρ¯h̃ctdV+ρ¯ũjhc̃njdS=(μ¯Prh̃cxjτhc,jsgs)njdS,ρ¯Ns̃tdV+ρ¯ũjÑsnjdS=(μ¯ScsNs̃xjτNs,jsgs)njdS+ṠNsdV,ρ¯Ys̃tdV+ρ¯ũjỸsnjdS=(μ¯ScsYs̃xjτYs,jsgs)njdS+ṠYsdV,
(6)

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

τϕ,jsgsρ¯(ϕuj̃ϕ̃ũj).
(7)

A gradient-diffusion closure model is used for τϕ,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, ṠZ, 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, hc, 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.

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

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,

sixiI(s)+(αa+αs)I(s)=αaσT4π+αs4πG.
(8)

Above, αa is the absorption coefficient, αs is the scattering coefficient, I(s) is the intensity along the direction si, T is the temperature, and the scalar flux is G. The black body radiation, Ib, is defined by σT4π. 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.

As noted in the introduction, the method of discrete ordinates45 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. 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,

Ωw(sixiI(s)+(αa+αs)IαaσT4παs4πG)dV+Ωw(sixiI(s)+(αa+αs)I)dV=0.
(9)

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:

Ωw(sixiI(s)+(αa+αs)IαaσT4παs4πG)dVΩIsiwxidV+ΓwsiInidS+Ωw(αa+αs)IdV=0.
(10)

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,

I=τh(sixiI(s)+(αa+αs)I(s)αaσT4παs4πG)=τhR(s).
(11)

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

Substituting Eq. (11) into Eq. (10) yields,

Ωw(sixiI(s)+(αa+αs)IαaσT4παs4πG)dV+ΩτhsiwxiR(s)dVΓτhwR(s)sinidSτhw(αa+αs)R(s)dV=0.
(12)

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

Ωw̃(sixiI(s)+(αa+αs)IαaσT4παs4πG)dVβΓτhwR(s)sinidS=0,
(13)

where w̃ is simply equal to,

w̃=w+τh(sjwxj+α(αa+αs)w).
(14)

When α=1, 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,

Ωw̃(xisiI(s)+(αa+αs)I(s)αaσT4παs4πG)dVβΓτhwR(s)sinidS=0,
(15)

and split into its Galerkin and stabilized contributions,

Ωw(xisiI(s)+(αa+αs)I(s)αaσT4παs4πG)dV+ΩτhsjwxjR(s)dV+αΩτhw(αa+αs)R(s)dVβΓτhwR(s)sinidS=0.
(16)

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

ΩwxisiI(s)dV=ΩI(s)siwxidV+ΓwsiI(s)nidS.
(17)

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:

Ω(I(s)siwxi+w[(αa+αs)I(s)αaσT4παs4πG])dV+ΓwsiI(s)nidS+ΩτhsjwxjR(s)dV+αΩτhw(αa+αs)R(s)dVβΓτhwR(s)sinidS=0.
(18)

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 = wI, 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,

wIxi=nIδ(xxΓI),
(19)

where ΓI is the boundary of control volume I and nI 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,

ΩI(s)sinidS+ΓsiI(s)nidS+Ω((αa+αs)I(s)αaσT4παs4πG)dVΩτhR(s)sinidS+αΩτh(αa+αs)R(s)dVβΓτhR(s)sinidS=0.
(20)

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 Burns56 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.,

I(s)=1π(τσT4+εσTw4+(1ετ)Hr).
(21)

Above, τ and ε are the transmissivity and emissivity, respectively, while the background reference temperature and surface temperature are T4 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

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 τadvh=h2 and τdiffh=1(αa+α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,

τh=(1(2h)2+(αa+αs)2)12.
(22)

In the above equation, the element length scale, h is determined by the dot product of the ordinate direction si 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 methodology56 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,

λTxjnj=ε(σT4G),
(23)

where λ is the thermal conductivity, nj 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.

FIG. 1.

(a) Concentric sphere geometric configuration and (b) code verification order-of-accuracy (L2 temperature norms).

FIG. 1.

(a) Concentric sphere geometric configuration and (b) code verification order-of-accuracy (L2 temperature norms).

Close modal

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 L2 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 (T10) was used while for the quadratic CVFEM basis, a quadrature order of 22 (T22) was activated. Sub-optimal convergence due to quadrature error pollution was found when using the T10 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 L2 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.

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.

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, χ=2DZ|Z|2 where DZ 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 χst=103s1) up to χst=50s1 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,

Yit=χ22YiZ2+ωiρ,
(24)
Tt=1ρcpi=1nsωihi+χ2(2TZ2+1cpTZcpZ+TZi=1nscp,icpYiZ)Hχmaxρcp1(1Zst)ZstTTTmax(t)T,
(25)

wherein Yi, hi, cp,i, and ωi are the mass fraction, enthalpy, isobaric specific heat capacity, and net mass production rate of species i, respectively, T is the temperature, t is the flamelet evolution time, ρ is the mass density, cp is the isobaric specific heat capacity of the mixture, T is a chosen cooled reference temperature, and 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 Pitch80 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 χst=0.1032s1. 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.

FIG. 2.

(a) Enthalpy and (b) temperature profiles for flamelets along the transient cooling trajectory for χst=0.1032s1. Lower profiles in each plot correspond to flamelet states later in evolution time.

FIG. 2.

(a) Enthalpy and (b) temperature profiles for flamelets along the transient cooling trajectory for χst=0.1032s1. Lower profiles in each plot correspond to flamelet states later in evolution time.

Close modal

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), (χ/2)2T/Z2, using an estimate that follows.81 

The change in the temperature gradient across the flame is proportional to (TmaxT)Zst1(1Zst)1, 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 (TmaxT), or equivalently the flame heat release, divided by a reaction zone thickness giving,

χ22TZ2|stχHZst(1Zst).
(26)

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 γ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 γst[0,2MJ/kg].

FIG. 3.

Maximum values of the temperature and OH mass fraction, and the stoichiometric enthalpy defect, over the dimensionless time, χstt, along cooling flamelet evolution trajectories for χst=0.0093,0.9627s1, with three distinct heat transfer coefficients, H=106,5×106,107W/m3.

FIG. 3.

Maximum values of the temperature and OH mass fraction, and the stoichiometric enthalpy defect, over the dimensionless time, χstt, along cooling flamelet evolution trajectories for χst=0.0093,0.9627s1, with three distinct heat transfer coefficients, H=106,5×106,107W/m3.

Close modal

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)T 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)T 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 γst is beginning to evolve toward a positive slope for the smallest 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.

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 Williams85) where molecular decomposition inhibits soot growth at high temperatures. To mimic this behavior, we introduce one additional term, ω̇R, that smoothly adjusts the soot nucleation and surface growth to a negative value above temperatures of Tc; the rate of this transition is given by an effective decomposition activation temperature, TR.

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 ω̇R, 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

ṠNs=ṠnucleationṠcoagulationYs̃1/6Ñs11/6,
(27)
ṠYs=WsṠnucleation+[ṠsurfacegrowthṠO2oxid.ṠOHoxid.]Ys̃2/3Ñs1/3,
(28)
ṠZ=ṠYs,
(29)
Ṡh=4σ(qsr¯Ys+qgasr¯)(αsoot¯Ys+αgas¯)G¯,
(30)

where αgas is determined from the mixture-averaged absorption coefficients of CO,CO2,CH4, and H2O, curve-fit from RADCAL data.86,87 Closures for soot production terms, e.g., ṠcoagulationYs̃1/6Ñs11/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, αsoot=max(α0+kαT,0)ρ/ρs with α0=375000m1,kα=1735m1K1, and ρs=1800 kg m−3.88 The radiation sources are qsr=αsootT4 and qgasr=αgasT4. The soot-related rates follow82 with the addition of the new ω̇R term

ω̇R=1exp(TR(1Tc1T)),
(31)
Ṡnucleation=AN[C2H2]exp(¯TN/T)ω̇R¯,
(32)
Ṡcoagulation=(24vdR/(ρsNA)))1/2ρ2T1/2¯,
(33)
Ṡsurfacegrowth=ASGvd2πρ[C2H2]exp(TSG/T)ω̇R¯,
(34)
ṠO2oxid.=AO2vd2πρ[O2]T1/2exp(TO2/T)¯,
(35)
ṠOHoxid.=AOHvd2πρ[OH]T1/2¯,
(36)

where TR=16000K,Tc=2050K,AN=54s1,TN=21100K,ASG=12000s1,TSG=12 000K,AO2=500s1,TO2=20000K,AOH=4.2325s1, and ρs=1800kgm3. Above, vd=61/3(πρsNA)1/3 converts volume to an effective particle diameter together with the Ys̃2/3Ñs1/3 to give a surface area dependence for the soot surface growth and oxidation terms; NA 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 Ws = 144 kg kmol−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.

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., Peters90) including ρ¯ 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

ϕ¯(Z̃,Z2̃,χ̃,γ̃)=001ϕ(Z,χ,γ)PZχγ(Z,χ,γ|Z̃,Z2̃,χ̃,γ̃)dZdχdγ,
(37)

where PZχγ is the joint PDF of the sub-filter fluctuations of ϕ in the Zχγ 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,

ϕ¯(Z̃,Z2̃,χ̃,γ̃)=01ϕ(Z,χ̃,γ̃)PZ(Z|Z̃,Z2̃)dZ,
(38)

in which PZ 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 ϕ̃ 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 ϕ̃ 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̃,Z2̃,χ̃, and γ̃, 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 χ̃,γ̃, and Z2̃ input values are converted to reference or scaled values χst,γst,Z2s, 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 η̃ represent the vector of filtered inputs acquired directly from evolved quantities in the simulation, and let η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, η̃table(ηo), and solving a minimization problem to identify the structured inputs that most closely reconstruct the filtered values:

ξ=argminη̃table(ξ)η̃,ϕ̃(η̃)=ϕ̃table(ξ).
(39)

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, χ(Z), the inversion of χ̃ to compute the corresponding stoichiometric value, χst, can be done explicitly. Similarly, an explicit conversion exists for the scaled scalar variance, Zs2=Z2̃/Z̃(1Z̃).

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 γst coordinate,

γst=argminγ̃table(γst)γ̃,ϕ̃(Z̃,Z2̃,χ̃,γ̃)=ϕ̃table(Z,Zs2,χst,γst),
(40)

where the Zs2 and χst values have been converted explicitly from the filtered input variables. We solve this equation for consistent values of γ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 γ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.

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 Shaddix93 due to the similarities between the mid- and large-scale pool fires.

The experimental configuration is modeled as a cylindrical domain 8D in diameter and 10D 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 1D to 1.5D) 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 10D 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).

FIG. 4.

Coarse mesh domain description outlining the (a) top-looking-down pool perspective and (b) the side view.

FIG. 4.

Coarse mesh domain description outlining the (a) top-looking-down pool perspective and (b) the side view.

Close modal

For turbulence modeling, the one-equation ksgs 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 (T4 = 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 type94 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 time95 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.

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.

FIG. 5.

Volume rendered temperature field for the (a) 5, (b) 2.5, and (c) 1.25 cm lip resolution.

FIG. 5.

Volume rendered temperature field for the (a) 5, (b) 2.5, and (c) 1.25 cm lip resolution.

Close modal

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ρ2ρ×p), 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.

FIG. 6.

Puffing cycle captured (roughly 0.143 s separation) in the order (a)–(l). Here, the reader can view the primary large-scale vortical structure that advects vertically in time and out of the image domain. This phenomena repeats itself.

FIG. 6.

Puffing cycle captured (roughly 0.143 s separation) in the order (a)–(l). Here, the reader can view the primary large-scale vortical structure that advects vertically in time and out of the image domain. This phenomena repeats itself.

Close modal

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,971.5D12. We note that uncertainty in the correlation is not provided although the authors report that “the agreement of the data is truly remarkable.”

FIG. 7.

Time signal post processing at a point two-diameters in height for (a) temperature and (b) the mixture fraction PSD. The primary peak noted at 0.60 Hz as compared to the correlation-based value of 0.67 Hz.

FIG. 7.

Time signal post processing at a point two-diameters in height for (a) temperature and (b) the mixture fraction PSD. The primary peak noted at 0.60 Hz as compared to the correlation-based value of 0.67 Hz.

Close modal

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 correlations41 is also shown. The correlation-based approach, which requires specification, predicts a plume height of 3.44D as compared to the simulation, 3.08D. 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.8D. 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.

FIG. 8.

(a) Centerline temperature and (b) velocity magnitude profiles. The Heskestad41 correlation is used to approximate the peak temperature location.

FIG. 8.

(a) Centerline temperature and (b) velocity magnitude profiles. The Heskestad41 correlation is used to approximate the peak temperature location.

Close modal

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

FIG. 9.

Surface radiation flux radial profile for three mesh resolutions; also included are experimental data and the first-generation RANS results.92 

FIG. 9.

Surface radiation flux radial profile for three mesh resolutions; also included are experimental data and the first-generation RANS results.92 

Close modal

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.

FIG. 10.

Temperature at varying normalized heights with different levels of refinement for (a) h/D = 0.5, (b) h/D = 1.0, (c) h/D = 2.0, and (d) h/D = 4.0.

FIG. 10.

Temperature at varying normalized heights with different levels of refinement for (a) h/D = 0.5, (b) h/D = 1.0, (c) h/D = 2.0, and (d) h/D = 4.0.

Close modal

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

FIG. 11.

Soot mass fraction at varying normalized heights with different levels of refinement for (a) h/D = 0.5, (b) h/D = 1.0, (c) h/D = 2.0, and (d) h/D = 4.0.

FIG. 11.

Soot mass fraction at varying normalized heights with different levels of refinement for (a) h/D = 0.5, (b) h/D = 1.0, (c) h/D = 2.0, and (d) h/D = 4.0.

Close modal

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, wk, the ordinate directions, sjk, and the local intensity, Ik, and is given by the expression, kwksjkIk. 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.

FIG. 12.

Magnitude of radiative heat flux at varying normalized heights with different levels of refinement for (a) h/D = 0.5, (b) h/D = 1.0, (c) h/D = 2.0, and (d) h/D = 4.0.

FIG. 12.

Magnitude of radiative heat flux at varying normalized heights with different levels of refinement for (a) h/D = 0.5, (b) h/D = 1.0, (c) h/D = 2.0, and (d) h/D = 4.0.

Close modal

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 Ahmed97 (0.6 Hz vs nominal 0.67 Hz). Peak temperature centerline values were predicted at a plume height of 3.08D 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.

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.

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.

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

1.
J.
Suo-Anttila
and
L.
Gritzo
, “
The effects of wind on fire environments containing large cylinders
,”
Combust. Sci. Technol.
181
,
68
77
(
2008
).
2.
L.
Gritzo
,
V.
Nicolette
,
S.
Tieszen
,
J.
Moya
, and
J.
Holen
, “
Heat transfer to the fuel surface in large pool fires
,” in
Transport Phenomena in Combustion
, edited by
S.
Chan
(
Taylor and Francis
,
Washington, DC
,
1995
).
3.
B.
Bainbridge
and
N.
Keltner
, “
Heat transfer to large objects in large pool fires
,”
J. Hazard. Mater.
20
,
21
40
(
1988
).
4.
S.
Tieszen
,
V.
Nicolette
,
L.
Gritzo
,
J.
Holen
,
D.
Murray
, and
J.
Moya
, “
Vortical structures in pool fires: Observation, speculation, and simulation
,”
Report No. SAND96-2607
(
Sandia National Laboratories SAND Series
,
Albuquerque, NM
,
1996
).
5.
S. R.
Tieszen
, “
On the fluids mechanics of fires
,”
Annu. Rev. Fluid Mech.
33
,
67
92
(
2001
).
6.
C.
Moen
,
S.
Tieszen
,
G.
Evans
,
A.
Black
,
S.
Domino
,
B.
Cochran
,
D.
Glaze
,
G.
Wagner
,
R.
Knaus
,
T.
Voskuilen
, and
F.
Pierce
, “
Sierra low-Mach module Fuego: Theory manual version 4.46
,”
Report No. SAND2017-10407
(
Sandia National Laboratories SAND Series
,
Albuquerque, NM
,
2017
).
7.
F.
Tamanini
, “
Reaction rates, air entrainment and radiation in turbulent fire plumes
,”
Combust. Flame
30
,
85
101
(
1977
).
8.
N.
Crauford
,
S.
Liew
, and
J.
Moss
, “
Experimental and numerical simulation of a bouyant fire
,”
Combust. Flame
61
,
63
77
(
1985
).
9.
V.
Novozhilov
and
H.
Koseki
, “
CFD predictions of pool fire burning rates and flame feedback
,”
Combust. Sci. Tech.
176
,
1283
1307
(
2004
).
10.
S.
Tieszen
,
S. P.
Domino
, and
A.
Black
, “
Validation of a simple turbulence model suitable for closure of temporally-filtered Navier-Stokes equations using a helium plume
,”
Report No. SAND2005-3210
(
Sandia National Laboratories SAND Series
,
Albuquerque, NM
,
2005
).
11.
S.
Tieszen
,
T.
O'Hern
,
E.
Weckman
, and
R.
Schefer
, “
Experimental study of the effect of fuel mass flux on a 1-m-diameter methane fire and comparison with a hydrogen fire
,”
Combust. Flame
139
,
126
141
(
2004
).
12.
A.
Luketa
,
V.
Romero
,
S.
Domino
,
D.
Glaze
,
M.
Sherman
, and
V.
Figueroa
, “
Validation and uncertainty quantification of Fuego simulations of calorimeter heating in a wind-driven hydrocarbon pool fire
,”
Report No. SAND2009-7605
(
Sandia National Laboratories SAND Series
,
Albuquerque, NM
,
2009
).
13.
A.
Brown
,
W.
Gill
,
G.
Evans
, and
D.
Jarboe
, “
Validation predictions of a 13 m/s cross-wind fire for Fuego and the University of Waterloo dataset
,”
Report No. SAND2008-0919
(
Sandia National Laboratories SAND Series
,
Albuquerque, NM
,
2008
).
14.
T.
Blanchat
,
P.
Helmick
,
R.
Jensen
,
A.
Luketa
,
R.
Deola
,
J.
Suo-Anttila
,
J.
Mercier
,
T.
Miller
,
A.
Ricks
,
R.
Simpson
,
B.
Demosthenous
,
S.
Tieszen
, and
M.
Hightower
, “
The Phoenix series large scale LNG pool fire experiments
,”
Report No. SAND2010-8676
(
Sandia National Laboratories SAND Series
,
Albuquerque, NM
,
2010
).
15.
L.
Gritzo
and
V.
Nicolette
, “
Coupling of large fire phenomenon with object geometry and object thermal response
,”
J. Fire Sci.
15
,
427
442
(
1997
).
16.
J.
Suo-Antilla
and
L.
Gritzo
, “
Thermal measurements from a series of tests with a large cylindrical calorimeter on the leeward edge of a JP-8 pool fire in crossflow
,”
Report No. SAND-20011986
(
Sandia National Laboratories
,
Albuquerque, NM
,
2001
).
17.
Y.
Wang
,
P.
Chatterjee
, and
J. L.
de Ris
, “
Large eddy simulation of fire plumes
,”
Proc. Combust. Inst.
33
,
2473
2480
(
2011
).
18.
H.
Wang
and
G.
Wang
, “
Interaction between crosswind and aviation-fuel fire engulfing a full-scale composite-type aircraft: A numerical study
,”
Aerospace
2
,
279
311
(
2015
).
19.
K.
McGrattan
,
H.
Baum
,
R.
Rehm
,
A.
Hamins
, and
G.
Forney
, “
Fire dynamics simulator technical reference guide
,”
Report No. NISTIR6467
(
National Institute of Standards and Technology Series
,
2007
).
20.
C.
Sedano
,
O.
Lopes
,
A.
Ladino
, and
F.
Munoz
, “
Prediction of a small-scale pool fire with FireFOAM
,”
Int. J. Chem. Eng.
2017
,
4934956
.
21.
G.
Maragkos
,
T.
Beji
, and
B.
Merci
, “
Advances in modelling in CFD simulations of turbulent gaseous pool fires
,”
Combust. Flame
181
,
22
38
(
2017
).
22.
P.
Chatterjee
,
J. L.
de Ris
,
Y.
Wang
, and
A. B.
Dorofeev
, “
A model for soot radiation in buoyant diffusion flames
,”
Proc. Combust. Inst.
33
,
2665
2671
(
2011
).
23.
P.
Chatterjee
,
Y.
Wang
,
K. V.
Meredith
, and
A. B.
Dorofeev
, “
Application of a subgrid soot-radiation model in the numerical simulation of a heptane pool fire
,”
Proc. Combust. Inst.
35
,
2573
2580
(
2015
).
24.
V. M.
Le
,
A.
Marchand
,
S.
Verma
,
R.
Xu
,
J.
White
,
A.
Marshall
,
T.
Rogaume
,
F.
Richard
,
J.
Luche
, and
A.
Trouve
, “
Simulations of a turbulent line fire with a steady flamelet combustion model coupled with models for non-local and local gas radiation effects
,”
Fire Saf. J.
106
,
105
113
(
2019
).
25.
S.
Domino
,
P.
Sakievich
, and
M.
Barone
, “
An assessment of atypical mesh topologies for low-Mach large-eddy simulation
,”
Comput. Fluids
179
,
655
669
(
2019
).
26.
L.
Jofre
,
S.
Domino
, and
G.
Iaccarino
, “
Eigensensitivity analysis of subgrid-scale stresses in large-eddy simulation of a turbulent axisymmetric jet
,”
Int. J. Heat Fluid Flow
77
,
314
335
(
2019
).
27.
L.
Jofre
,
S.
Domino
, and
G.
Iaccarino
, “
A framework for characterizing structural uncertainty in large-eddy simulation closures
,”
Flow Turbul. Combust.
99
,
1
23
(
2017
).
28.
S.
Domino
, “
Design-order, non-conformal low-Mach fluid algorithms using a hybrid CVFEM/DG approach
,”
J. Comput. Phys.
359
,
331
351
(
2018
).
29.
P.
Lin
,
M.
Bettencourt
,
S.
Domino
,
T.
Fisher
,
M.
Hoemmen
,
J.
Hu
,
E.
Phipps
,
A.
Prokopenko
,
S.
Rajamanickam
,
C.
Siefert
, and
S.
Kennon
, “
Towards extreme-scale simulations for low-Mach fluids with second-generation Trilinos
,”
Parallel Proccess. Lett.
24
,
1442005
(
2014
).
30.
S.
Domino
and
M.
Barone
, “
Wind turbine blade wall-resolved large-eddy simulation
,”
Report No. SAND2016-4085C
(
Sandia National Laboratories SAND Series
,
Albuquerque, NM
,
2016
).
31.
E.
Weckman
and
A.
Strong
, “
Experimental investigation of the turbulence structure of medium-scale methanol pool fire
,”
Combust. Flame
105
,
245
266
(
1996
).
32.
V.
Nicolette
,
S.
Tieszen
,
A.
Black
,
S. P.
Domino
, and
T.
O'Hern
, “
A turbulence model for buoyant flows based on vorticity generation
,”
Report No. SAND2005-6273
(
Sandia National Laboratories SAND Series
,
Albuquerque, NM
,
2005
).
33.
M. M.
Ahmed
and
A.
Trouvé
, “
Large eddy simulation of the unstable flame structure and gas-to-liquid thermal feedback in a medium-scale methanol pool fire
,”
Combust. Flame
225
,
237
254
(
2021
).
34.
L.
Ma
,
F.
Nmira
, and
J.-L.
Consalvi
, “
Large eddy simulation of medium-scale methanol pool fires—Effects of pool boundary conditions
,”
Combust. Flame
222
,
336
354
(
2020
).
35.
H.
Koo
,
J.
Hewson
, and
R.
Knaus
, “
LES soot-radiation predictions of buoyant fire plumes
,”
Report No. SAND2018-3006C
(
Sandia National Laboratories SAND Series
,
Albuquerque, NM
,
2018
).
36.
W.
Yao
,
J.
Zhang
,
A.
Nadjai
,
T.
Beji
, and
M. A.
Delichatsios
, “
A global soot model developed for fires: Validation in laminar flames and application in turbulent pool fires
,”
Fire Saf. J.
46
,
371
387
(
2011
).
37.
Z.
Chen
,
J.
Wen
,
B.
Xu
, and
S.
Dembele
, “
Extension of the eddy dissipation concept and smoke point soot model to the les frame for fire simulations
,”
Fire Saf. J.
64
,
12
26
(
2014
).
38.
See http://uintah.utah.edu/projects/15-projects/projects-old/16-csafe.html for “
Center for Simulation of Accidental Fires and Explosions
(
2020
)” (last accessed February, 2021).
39.
J.
Spinti
,
J.
Thornock
,
E.
Eddings
,
P.
Smith
, and
A.
Sarofim
,
Transport Phenomena in Fires: Heat Transfer to Objects in Pool Fires
, WIT Transactions on State-of-the-Art in Science and Engineering (
WIT
,
Billerica, MA
,
2008
), pp.
69
136
.
40.
T.
O'Hern
,
E.
Weckman
,
A.
Gerhart
,
S.
Tiezen
, and
R.
Sheffer
, “
Experimental study of a turbulent buoyant helium plume
,”
J. Fluid Mech.
544
,
143
171
(
2005
).
41.
G.
Heskestad
, “
Fire plumes
,” in
SFPE Handbook of Fire Protection Engineering
, edited by
P.
DiNenno
(
1995
), pp.
2:9
2:19
.
42.
J.
de Ris
, “
Mechanism of buoyant turbulent diffusion flames
,”
Procedia Eng.
62
,
13
27
(
2013
).
43.
J.
Tencer
and
J.
Howell
, “
Coupling radiative heat transfer in participating media with other heat transfer modes
,”
J. Braz. Soc. Mech. Sci. Eng.
38
,
1473
1487
(
2016
).
44.
K.
Jensen
,
J.
Ripoll
,
D.
Wray
,
A.
Joseph
, and
M.
El Hafi
, “
Radiative transfer modeling of a large pool fire by discrete ordinates, discrete transfer, ray tracing, Monte Carlo, and moment methods
,” in
Studying Turbulence Using Numerical Simulation Databases - X
, edited by
P.
Moin
and
N.
Mansour
(
Stanford Center for Turbulence Research
,
2004
), pp.
353
366
.
45.
M.
Modest
,
Radiative Heat Transfer
, 3rd ed. (
Academic Press
,
Oxford, UK
,
2013
).
46.
P. J.
Coelho
, “
Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media
,”
J. Quant. Spectrosc. Radiat. Transfer
145
,
121
146
(
2014
).
47.
J.
Chai
,
H.
Lee
, and
S.
Patankar
, “
Finite volume method for radiation heat transfer
,”
J. Thermophys. Heat Transfer
8
,
419
425
(
1994
).
48.
J.
Strohle
,
U.
Schnell
, and
K.
Hein
, “
A mean flux discrete ordinates interpolation scheme for general coordinates
,” in
3rd International Conference on Heat Transfer
(
AIAA
,
Washington, DC
,
2001
).
49.
M.
Sakami
,
A.
Charette
, and
V.
LeDez
, “
Radiative heat transfer in 3-dimensional enclosures of complex geometry using the discrete ordinates method
,”
J. Quant. Spectrosc. Radiat. Transfer
59
,
117
136
(
1998
).
50.
J.
Liu
,
H.
Shang
,
Y.
Chen
, and
T.
Wang
, “
Development of an unstructured radiation model applicable for two dimensional planar, axisymmetric and 3-dimensional geometries
,”
J. Quant. Spectrosc. Radiat. Transfer
66
,
17
33
(
2000
).
51.
D.
Joseph
,
M.
El Hafi
,
R.
Fournier
, and
B.
Cuenot
, “
Comparison of three spatial differencing schemes in discrete ordinates method using three-dimensional unstructured meshes
,”
Int. J. Therm. Sci.
44
,
851
864
(
2005
).
52.
W.
Fiveland
and
J.
Jessee
, “
Comparison of discrete ordinates formulations for radiative heat transfer in multidimensional geometries
,”
J. Thermophys. Heat Transfer
9
,
47
54
(
1995
).
53.
R.
Castro
and
J.
Trelles
, “
Spatial and angular finite element method for radiative transfer in participating media
,”
J. Quant. Spectrosc. Radiat. Transfer
157
,
81
105
(
2015
).
54.
W.
Fiveland
and
J.
Jessee
, “
Finite element formulation of the discrete-ordinates method for multidimensional geometries
,”
J. Thermophys. Heat Transfer
8
,
171
180
(
1994
).
55.
M.
Avila
,
R.
Codina
, and
J.
Principe
, “
Spatial approximation of the radiation transport equation using a two-scale finite element method
,”
Comput. Methods Appl. Mech. Eng.
200
,
425
438
(
2011
).
56.
S.
Burns
,
Application of Spatial and Angular Domain Based Parallelism to a Discrete Ordinates Formulation with Unstructured Spatial Discretization
(
International Centre for Heat and Mass Transfer
,
1997
).
57.
H.
Edwards
,
A.
Williams
,
G.
Sjaardema
,
D.
Baur
, and
W.
Cochran
, “
Sierra toolkit computational mesh computational model
,”
Report No. SAND-20101192
(
Sandia National Laboratories
,
Albuquerque, NM
,
2010
).
58.
M.
Heroux
,
R.
Bartlett
,
V.
Howle
,
R.
Hoekstra
,
J.
Hu
,
T.
Kolda
,
R.
Lehoucq
,
K.
Long
,
R.
Pawlowski
,
E.
Phipps
,
A.
Salinger
,
J.
Thornquist
,
R.
Tuminaro
,
J.
Willenbring
, and
A.
Williams
, “
An overview of Trilinos
,”
Report No. SAND-20032927
(
Sandia National Laboratories
,
Albuquerque, NM
,
2003
).
59.
S.
Domino
, “
Sierra low-Mach module Nalu: Theory manual
,”
Report No. SAND2015-3107W
(
Sandia National Laboratories SAND Series
,
Albuquerque, NM
,
2014
).
60.
S.
Domino
,
G.
Wagner
,
A.
Luketa-Hanlin
,
A.
Black
, and
J.
Sutherland
, “
Verification for multi-mechanics applications
,” AIAA Paper No. 2007-1933,
2007
.
61.
See https://www.energy.gov for “
Liquefied natural gas safety research: Report to Congress
(
2012
)” (last accessed December, 2020).
62.
S.
Rodrigous
,
S.
Domino
, and
M.
El-Genk
, “
Fluid flow and heat transfer an analysis of the VHTR lower plenum using the Fuego CFD code
,”
Computational Fluid Dynamics (CFD) for Nuclear Reactor Safety Applications (CFD4NRS-3)
(
2010
).
63.
S.
Domino
,
F.
Pierce
, and
J.
Hubbard
, “
A multi-physics computational investigation of droplet pathogen transport emanating from synthetic coughs and breathing
,”
Atomization Spray
31
,
1
23
(
2021
).
64.
S. P.
Domino
, “
A case study on pathogen transport, deposition, evaporation and transmission: Linking high-fidelity computational fluid dynamics simulations to probability of infection
,”
Int. J. Comput. Fluid Dyn.
0
,
1
15
(
2021
).
65.
S.
Paolucci
, “
On the filtering of sound waves from the Navier-Stokes equations
,”
Report No. SAND82-8257
(
Sandia National Laboratories SAND Series
,
Albuquerque, NM
,
1982
).
66.
M.
Germano
,
U.
Piomelli
,
P.
Moin
, and
W. H.
Cabot
, “
A dynamic subgrid-scale eddy viscosity model
,”
Phys. Fluids
3
,
1760
1765
(
1991
).
67.
A.
Yoshizawa
, “
Bridging between eddy-viscosity-type and second-order models using a two-scale DIA
,”
Int. Symp. Turbl. Shear Flow
3
,
23.1.1
23.1.6
(
1993
).
68.
S.
Domino
, “
Toward verification of formal time accuracy for a family of approximate projection methods using the method of manufactured solutions
,” in
Studying Turbulence Using Numerical Simulation Databases - XI
, edited by
P.
Moin
and
N.
Mansour
(
Stanford Center for Turbulence Research
,
2006
), pp.
163
177
.
69.
M. J.
Zimberg
,
S. H.
Frankel
,
J. P.
Gore
, and
Y. R.
Sivathanu
, “
A study of coupled turbulent mixing, soot chemistry, and radiation effects using the linear eddy model
,”
Combust. Flame
113
,
454
469
(
1998
).
70.
A. J.
Ricks
,
J. C.
Hewson
,
A. R.
Kerstein
,
J. P.
Gore
,
S. R.
Tieszen
, and
W. T.
Ashurst
, “
A spatially developing one-dimensional turbulence (ODT) study of soot and enthalpy evolution in meter-scale buoyant turbulent flames
,”
Combust. Sci. Technol.
182
,
60
101
(
2010
).
71.
G.
Schneider
and
M.
Raw
, “
A skewed, positive influence coefficient upwinding procedure for control-volume-based finite-element convection-diffusion computation
,”
Numer. Heat Transfer
9
,
1
26
(
1986
).
72.
S.
Domino
, “
A comparison between low-order and higher-order low-Mach discretization approaches
,” in
Studying Turbulence Using Numerical Simulation Databases - XV
, edited by
M.
P.
and
J.
Urzay
(
Stanford Center for Turbulence Research
,
2014
), pp.
387
396
.
73.
C.
Thurgood
,
A.
Pollard
, and
A.
Becker
, “
The TN quadrature set for the discrete ordinates method
,”
J. Heat Transfer
117
,
1068
1070
(
1995
).
74.
T.
Hughes
,
G.
Scovazzi
, and
L.
Franca
, “
Multiscale and stabilized methods
,” in
Encyclopedia of Computational Mechanics
, edited by
E.
Stein
,
R.
de Borst
, and
T.
Hughes
(
John Wiley and Sons
,
2004
), Chap. 2, pp.
5
59
.
75.
M.
Martinez
, “
Comparison of Galerkin and control volume finite element for advection-diffusion problems
,”
Int. J. Numer. Methods Fluids
50
,
347
376
(
2006
).
76.
C. R.
Swaminathan
and
V. R.
Voller
, “
Streamwise upwind scheme for control-volume finite elements, Part I. Formulations
,”
Numer. Heat Transfer B
22
,
95
107
(
1992
).
77.
G.
Blanquart
,
P.
Pepiot-Desjardins
, and
H.
Pitsch
, “
Chemical mechanism for high temperature combustion of engine relevant fuels with emphasis on soot precursors
,”
Combust. Flame
156
,
588
607
(
2009
).
78.
H.
Pitsch
,
M.
Chen
, and
N.
Peters
, “
Unsteady flamelet modeling of turbulent hydrogen-air diffusion flames
,”
Symp. (Int.) Combust.
27
,
1057
1064
(
1998
).
79.
H.
Pitsch
, “
Unsteady flamelet modeling of differential diffusion in turbulent jet diffusion flames
,”
Combust. Flame
123
,
358
374
(
2000
).
80.
M.
Ihme
and
H.
Pitsch
, “
Modeling of radiation and nitric oxide formation in turbulent nonpremixed flames using a flamelet/progress variable formulation
,”
Phys. Fluids
20
,
055110
(
2008
).
81.
J. C.
Hewson
, “
An extinction criterion for nonpremixed flames subject to brief periods of high dissipation rates
,”
Combust. Flame
160
,
887
897
(
2013
).
82.
I.
Aksit
and
J.
Moss
, “
A hybrid scalar model for sooting turbulent flames
,”
Combust. Flame
145
,
231
244
(
2006
).
83.
F.
Bisetti
,
G.
Blanquart
,
M. E.
Mueller
, and
H.
Pitsch
, “
On the formation and early evolution of soot in turbulent nonpremixed flames
,”
Combust. Flame
159
,
317
335
(
2012
).
84.
M.
Frenklach
,
S.
Taki
,
M. B.
Durgaprasad
, and
R. A.
Matula
, “
Soot formation in shock tube pyrolysis of acetylene, allene and 1,3-butadiene
,”
Combust. Flame
54
,
81
101
(
1983
).
85.
A.
Alexiou
and
A.
Williams
, “
Soot formation in shock-tube pyrolysis of toluene-n-heptane and toluene-iso-octane mixture
,”
Fuel
74
,
153
158
(
1995
).
86.
W. L.
Grosshandler
, “
‘Radcal: A narrow-band model for radiation,’ Calculations in a Combustion Environment
,”
Report No. 1402
(NIST,
1993
).
87.
R.
Barlow
,
A.
Karpetis
,
J.
Frank
, and
J.-Y.
Chen
, “
Scalar profiles and no formation in laminar opposed-flow partially premixed methane/air flames
,”
Combust. Flame
127
,
2102
2118
(
2001
).
88.
G. L.
Hubbard
and
C. L.
Tien
, “
Infrared mean absorption coefficient of luminous flames and smoke
,”
J. Heat Transfer
100
,
235
239
(
1978
).
89.
K. M.
Leung
,
R. P.
Lindstedt
, and
W. P.
Jones
, “
A simplified reaction mechanism for soot formation in nonpremixed flames
,”
Combust. Flame
87
,
289
305
(
1991
).
90.
N.
Peters
,
Turbulent Combustion
(
IOP Publishing
,
Cambridge, UK
,
2001
).
91.
E.
Eddings
,
Y.
Shihong
,
C.
William
, and
A.
Sarafim
, “
Formulation of a surrogate for the simulation of jet fuel pool fires
,”
Combust. Sci. Technol.
177
,
715
739
(
2005
).
92.
S.
Domino
,
C.
Moen
,
S.
Burns
, and
G.
Evans
,
SIERRA/Fuego: A Multi-Mechanics Fire Environment Simulation Tool
(
AIAA
,
Washington, DC
,
2003
).
93.
J.
Murphy
and
C.
Shaddix
, “
Soot property measurements in a two-meter diameter JP-8 pool fire
,”
Combust. Sci. Technol.
178
,
865
894
(
2006
).
94.
D. S.
Burgess
,
A.
Strasser
, and
J.
Grumer
,
Diffusive Burning of Liquid Fuels in Open Trays
(
American Chemical Society
,
Washington, DC
,
1961
).
95.
V.
Babrauskas
, “
Estimating large pool fire burning rates
,”
Fire Technol.
19
,
251
261
(
1983
).
96.
K.
Fukumoto
,
J. X.
Wen
,
M.
Li
,
Y.
Ding
, and
C.
Wang
, “
Numerical simulation of small pool fires incorporating liquid fuel motion
,”
Combust. Flame
213
,
441
454
(
2020
).
97.
T.
Cetegen
and
B.
Ahmed
, “
Experiments on the periodic instability of buoyant plumes and pool fires
,”
Combust. Flame
93
,
157
184
(
1993
).
98.
T.
Blanchat
,
L.
Humphries
, and
W.
Gill
, “
Sandia heat flux gauge thermal response and uncertainty models
,”
Report No. SAND2000-1111
(
Sandia National Laboratories SAND Series
,
Albuquerque, NM
,
2000
).
99.
T.
Blanchat
and
V.
Figueroa
, “
Large-scale open pool experimental data and analysis for fire model validation and development
,”
Fire Saf. Sci.
9
,
105
115
(
2008
).