Kirk A Jensen
email: jensenk@pixia.com.

Jean-François Ripoll
email: ripoll@lanl.gov.

Alan A Wray

David Joseph

Mouna El

On various modeling approaches to radiative heat transfer in pool fires

Keywords: Fire, Radiation, Discrete ordinates, Discrete transfer, Monte Carlo, Ray tracing, M 1 model

Six computational methods for solution of the radiative transfer equation in an absorbing-emitting, nonscattering gray medium were compared for a 2-m JP-8 pool fire. The emission temperature and absorption coefficient fields were taken from a synthetic fire due to the lack of a complete set of experimental data for computing radiation for large and fully turbulent fires. These quantities were generated by a code that has been shown to agree well with the limited quantity of relevant data in the literature. Reference solutions to the governing equation were determined using the Monte Carlo method and a ray-tracing scheme with high angular resolution. Solutions using the discrete transfer method (DTM), the discrete ordinates method (DOM) with both S 4 and LC 11 quadratures, and a moment model using the M 1 closure were compared to the reference solutions in both isotropic and anisotropic regions of the computational domain. Inside the fire, where radiation is isotropic, all methods gave comparable results with good accuracy. Predictions of DTM agreed well with the reference solutions, which is expected for a technique based on ray tracing. DOM LC 11 was shown to be more accurate than the commonly used S 4 quadrature scheme, especially in anisotropic regions of the fire domain. On the other hand, DOM S 4 gives an accurate source term and, in isotropic regions, correct fluxes. The M 1 results agreed well with other solution techniques and were comparable to DOM S 4 . This represents the first study where the M 1 method was applied to a combustion problem occurring in a complex three-dimensional geometry. Future applications of M 1 to fires and similar problems are recommended, considering its similar accuracy and the fact that it has significantly lower computational cost than DOM S 4 .

Introduction

Accurate prediction of the heat transfer from a large hydrocarbon fire, which can occur from an industrial or transportation accident, is important for consideration of the thermal hazard to engineered systems, personnel, and facilities. Fires of this scale typically have relatively low velocities and high temperatures and produce significant quantities of soot. The majority of heat transfer in such a fire is dominated by radiative emission from the high-temperature soot [1]. As a consequence, accurate solution of the radiative transfer equation is important for prediction of the resulting thermal hazard to an object exposed to the fire.

The cost and accuracy of a large fire simulation are strongly dependent on the choice of numerical model for solving the radiative transfer equation (RTE) governing this phenomenon. The RTE is complicated by the fact that, in addition to the three-dimensional space variables and time commonly found in fluids problems, integration over all directions of propagation is necessary at each point in the domain. The radiative source term, which couples to the energy equation, must be computed with sufficient accuracy to ensure a correct prediction of the evolution of the fire. Moreover, the scale of these fires, on the order of a few meters, as well as the existence of large spatial and temporal gradients in physical properties on microscales, would require simulations using millions/billions of grid points and small time steps, making it impractical to routinely solve the governing equations directly.

Approximations, sometimes drastic, in terms of models or coarse meshes used for the radiative transfer are required, as is often done for turbulence, chemistry, and soot production. For example, the work of Porterie and Loraud [2,3] demonstrates the number of assumptions and approximations needed for modeling the radiation in compartment fires. The most common method of accelerating the numerical solution of the radiant transport is to limit the number of rays of integration. However, the choice of the minimum number of angles is often not rigorously justified, and the risk of choosing a poor resolution exists in decreasing the computational time (see [4] and references therein). The challenge, therefore, is choosing a computationally efficient numerical solution method that predicts radiation in fires, as well as radiative flux to objects inside and/or outside the fire, with sufficient accuracy.

For comprehensive coverage of the physics of fire and current research in fire modeling, the authors refer the reader to the following recent reviews: Tieszen [5] and Drysdale [START_REF] Drysdale | An Introduction to Fire Dynamics[END_REF] for fire dynamics; Joulain [START_REF] Joulain | [END_REF] for pool fires; Novozhilov [8] concerning modeling and computation of compartment fires; and, more recently, Sacadura [4] on radiation in fires.

In this study, six common numerical methods for solving the radiative transfer equation were compared when applied to a realistic synthetic full-field, threedimensional fire (described below). The six methods include ray tracing, discrete transfer (DTM), discrete ordinates (DOM), moment (with the M 1 closure), and Monte Carlo (MC). Ray tracing and the discrete transfer method are based on straightforward and direct integration of the integral equation by tracing a specified number of rays originating from each point throughout the domain. In contrast, the DOM consists of solving transport equations by finite volume methods along discrete directions. The angular integration in DOM is performed with a select numerical quadrature scheme. In this study, two quadrature schemes were investigated, including the traditional S 4 scheme with 24 angles, and the recently investigated LC 11 formulation of Lathrop and Carlson with 96 angles [9]. The Monte Carlo method used here is formulated in terms of net exchange and presumes the form of the probability density function for efficient computations. The M 1 moment method in this study consists of a set of four hyperbolic equations obtained by integration of the RTE over frequency and angles. This method, however, requires a closure model for the radiative pressure. The maximum entropy closure [10] is used here. Our study is unfortunately limited to the only methods which were available to us.

General descriptions of these methods can be found in [START_REF] Modest | Radiative Heat Transfer[END_REF][START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF][START_REF] Mihalas | Foundation of Radiation Hydrodynamics[END_REF][START_REF] Pomraning | The Equations of Radiation Thermodynamics[END_REF], and more details will be covered below with comprehensive references describing the implementations. The discrete transfer and discrete ordinates methods are commonly used to solve radiation in fires. For example, the DTM is utilized in Vulcan [START_REF] Holen | Proc. Combust. Inst[END_REF][START_REF] Brown | ASME Summer Heat Transfer Conference[END_REF], JASMINE [START_REF] Markatos | [END_REF], and SOFIE [18], and DOM is used in various French fire codes [2,3,[START_REF] Joulain | Proc. Int. Symp. Fire Safety Sci[END_REF] and SIERRA/Syrinx [START_REF] Burns | [END_REF][START_REF] Moen | Proceeding of IMECE2002-33098[END_REF]. The Monte Carlo method has been used for computation of radiation in compartment fires [START_REF] Snegirev | [END_REF], although it is not generally used in fire codes, given its computational cost. The M 1 method has not been applied previously to a threedimensional fire problem to the authors' knowledge.

To compare the accuracy of the commonly used DTM and DOM techniques, as well as the M 1 method, the solution of the heat flux and radiative source terms are required. Unfortunately, a data set of these quantities is unavailable for this purpose, as explained in the next section. Therefore, the ray tracing and Monte Carlo techniques were used to obtain these reference-solution fields. With the "known" radiative fields obtained from Monte Carlo and ray tracing, the accuracy of the DT, DO, and M 1 methods can be quantified. Both the effects of different angular res-olutions and general accuracy of the computational scheme were addressed in this study.

This represents the first study in which the M 1 moment formulation is applied to the computation of fires. The closure approximation was first introduced by Minerbo [10] while searching for an appropriate Eddington factor to close the radiative pressure term, and has been developed over the past 20 years by the efforts of [23][24][25] and numerous works cited in [26]. However, it has received little attention in the combustion community and little theoretical description of M 1 exists in modern radiation textbooks [START_REF] Modest | Radiative Heat Transfer[END_REF][START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF]. Consequently, the implementation details for application to fires will be discussed below. This method is particularly attractive because it has the lowest theoretical computational cost compared to the other methods. The angular dependence in M 1 is handled analytically, thus removing two variables from the integration of the gray RTE (leaving time and position as the remaining variables).

The paper is constructed as follows. The motivation for and creation of the synthetic fire are described, including the geometry of the facility it represents, the computation of the fire, and the boundary conditions used. This description is followed by the definition of the radiant transport problem and the presentation of the six computation methods used to solve the radiation in this fire field. Discussion of the numerical results focuses on two aspects of solving the radiative transport problem of interest to fire modeling. First, an angular study was performed and the effects of choosing too low an angular resolution on the coupling to the fire hydrodynamics are considered with the ray tracing and discrete transfer methods. Second, the six methods were used to solve the radiative transfer of the fire, and particular attention was given to both the radiative flux and the radiative source term profiles. The reference solutions from Monte Carlo and ray tracing were used to quantify and compare the accuracy of these variables obtained by the discrete transfer, discrete ordinates, and M 1 moment methods.

Synthetic fire

Solution of radiant transport in a hydrocarbon fire requires knowledge of the absorption coefficient and emission temperature fields throughout the participating medium. Those coefficients involve knowledge of the soot, CO 2 , and H 2 O concentrations and the temperature. Ideally, the radiative properties would be fully spatially resolved so that no filtering of the radiative transport equation was required. Unfortunately, such a data set is nowhere near obtainable. For instance, experimental measurements as detailed as a direct numerical simulation (DNS) of the turbulent reacting flow field would require between 10 12 and 10 15 measurement points (corresponding to between 200 and 20 µm needed for measuring the soot volume fraction accurately within the gradient of temperature in a 2-m fire).

The goal of the current study is not to assess radiative transport methods for DNS solutions of large fires, for the same reason that data sets do not existthey are not obtainable. Current and foreseeable computational strategies for large fires will require that the radiative transport equations (along with the turbulent reacting flow equations) be modeled. The filtered radiative transport equation then contains effective emission temperature and absorption properties that account for the high spatial frequencies that are not resolvable on the grid. These models are sometimes referred to as turbulence-radiation interaction models (TRI). Typically, correlations of the form Y s T 5 , where Y s is the soot volume fraction and T the temperature, are needed as soon as "low resolution" profiles are used. (The reader should note here that using a fourth power of T in the emission term is not appropriate for large fires because it neglects TRI and leads to large error.) Hence fire temperatures, for which there are many data sets, cannot be used as a validation field.

Unfortunately, data sets for even filtered properties are not available either, except for limited measurements along the centerline and a radial [27], not enough for assessing three-dimensional radiation transport methods. The lack of measurements is due to complex fire dynamics [28] and difficulties associated with diagnostic instruments in high-temperatures sooting environments [1,29,30]. Given this lack of comprehensive experimental data, a synthetic 2-m JP-8 fire was created with the Vulcan fire simulation tool from Sandia Albuquerque [START_REF] Holen | Proc. Combust. Inst[END_REF]. Vulcan has been used in recent years to simulate pool fires resembling those in the Fire Laboratory for Accreditation of Models and Experiments (FLAME) facility in Albuquerque, New Mexico [START_REF] Brown | ASME Summer Heat Transfer Conference[END_REF].

FLAME consists of a 6.1-m square water-cooled wall structure with a 6-m air ring on the floor and a 2-m fuel pan elevated approximately 1.6 m above ground. Fig. 1 shows (left) the exterior of the FLAME facility and (right) a video still image during a 2-m JP-8 fire. The heights above the fuel pan are indicated by the horizontal bars in the figure, and the scale of the fuel pan is included at the bottom. The 3D computational mesh of this facility consists of a 92 × 92 × 120 Cartesian grid and its central plane contour plot is shown in Fig. 2 (left). The dark cells in this figure represent the walls, the air inlet ring (bottom), the fuel pan (center), and the atmosphere outside the chimney (top outer edges). It has been confirmed experimen- tally [START_REF] Blanchat | SAND Report[END_REF] that the four corners of the square facility do not influence the axial symmetry of the source air flow, and the resulting fires are indeed axially symmetric. The calculations were performed using the full 3D data set as input. The symmetry of all solutions was verified and, therefore, results are only presented for various axial and radial locations along the central plane of the facility.

To compute the fire, the fluid conservation equations are solved using a Reynolds-averaged Navier-Stokes (RANS) approach, with sufficient iterations from ignition to reach a steady-state solution. The emission temperature, defined as ( T 4 ) 1/4 , is plotted in Fig. 1 (right) and accounts for TRI. To obtain a relevant estimate of the emission temperature, a simple model for turbulent radiation interactions is used employing a double delta function representing soot in the flame sheets and soot that has been mixed into the cooler surroundings [START_REF] Holen | Proc. Combust. Inst[END_REF]. The delta function weights are related to an estimate of the volume fraction of flame in each cell.

Vulcan treats the fire as a gray absorbing-emitting medium with negligible scattering [START_REF] Houf | SAND Report[END_REF]. The absorption coefficients contain contributions mostly from soot particles but also include contributions from carbon dioxide and water vapor absorption at select wavelengths with an empirically based wideband model. The model for these coefficients assumes that the medium is gray since soot is the dominant absorbing and emitting species [START_REF] Kearney | Proc. ASME Int. Mechanical Engineering Congress and Exposition[END_REF] as we now explain. Soot is a broadband emitter the spectral absorption of which follows a linear spectral law on the whole frequency spectrum (see for instance [START_REF] Lee | [END_REF]). This law is evaluated here with the coefficients of [35] and is integrated over the whole frequency spectrum in order to obtain a Planck mean coefficient which is used afterward in the RTE. Thus, the nongray character of soot is accounted for through this Planck mean coefficient. The filtered absorption coefficient is plotted in Fig. 3.

Normally in the solution of the RTE, it is necessary to know the optical properties of all surfaces with which radiative energy can come in contact. At a surface, all incident energy must be absorbed, transmitted, or reflected-that is, α + τ + ρ = 1, where α, τ , and ρ are the surface absorptance, transmittance, and reflectance, respectively. When the domain of computation is as complex as the FLAME facility (see Fig. 2), implementing the prescribed boundary conditions uniformly across five different codes proves challenging. Instead, an alternative approach using "ghost cells," which can be considered similar to the immersed boundary conditions technique, was used to overcome this difficulty. The domain is considered as an open domain throughout which radiation can travel. The walls and the pan are simulated by cells having a constant and very large opacity (a few hundred m -1 ) and an ambient temperature (293 K), such that the walls and the pan are artificially maintained in a black and cold state, conditions that hold in the FLAME facility. At the exhaust opening at the top of the chimney, it is assumed that all energy leaves unimpeded (τ = 1).

Solution methods of radiative transfer

The gray radiative transfer equation (RTE) describes the change in radiation intensity, I , through an absorbing and emitting gray medium along a path of length ds in a solid angle Ω defined around the direction of propagation s [START_REF] Modest | Radiative Heat Transfer[END_REF], For most heat transfer applications, the primary engineering quantities of interest are the net incident radiation (G), the radiative flux (q r ), and the divergence of the radiative heat flux (∇ • q r ), called the radiative source term hereafter. These quantities can be derived from the following integrals of intensity over solid angle:

G = 4π I (s) dΩ, q r = 4π I (s)s dΩ, and (2)
∇ • q r = κ(4πI b -G).
Solution of the RTE (Eq. ( 1)) and Eqs. (2) using each solution scheme is outlined separately below.

Discrete transfer method (DTM)

The discrete transfer method (DTM) used in Vulcan is an enhanced version of the original model proposed by Shah [START_REF] Shah | The Computation of Radiation Heat Transfer[END_REF][START_REF] Lockwood | [END_REF]. The enhancements were selected to obtain an acceptable compromise between accuracy and calculation speed. The technique will be tested by comparing the results obtained herein with those obtained from verified and highly accurate Monte Carlo and ray-tracing techniques.

Within the computational domain a radiation box is defined to speed the calculation by focusing on the region with high thermal emission. This box defines where rays originate in the tracing technique. For this study, the box was defined as the smallest grid-conformal parallelepiped encompassing all control volumes with a emission temperature greater than 800 K.

For each node on the boundary of the box, a specified number of rays are emitted over a hemisphere and followed to the boundary of the calculation domain; a corresponding ray is followed back from the boundary to the original point, and on to the next boundary. Along these traces, the change of intensity from absorption and emission is calculated over each control volume in the path with proper weighting given to the solid angle and the originating projected area.

The change of intensity for the ray within a control volume is found from a recurrence relation obtained from analytical integration of Eq. ( 1), ( 3)

I n+1 = I n exp(-κδs) + I b 1 -exp(-κδs) ,
where δs is the distance over which the beam passes through the control volume.

The source term for the energy equations, Eqs. (2), is found by summing the net gain or loss of radiation energy in each control volume intersected during a ray trace. The contribution to the source term from one beam i passing through a control volume n is given by (4) S n,i = (I n+1 -I n )s i dA dΩ, where dA is the area from the element at the ray origin boundary and Ω i is the solid angle represented by the beam. The total radiant source term for the nth control volume is found by summing over N total beams:

(5) Q r dV = i=1,N S n,i .
The heat flux to a surface is not explicitly calculated throughout the field of Vulcan, but rather at selected surfaces (i.e., cell faces). The hemispherical flux in W/m 2 is derived from this model by integrating all incoming rays on a surface. This integration requires a large number of rays to be traced from each node of the radiation box to be accurate, but from experience it has been found to be quite fast when a limited number of selected surfaces are used. To compare to the other methods, the hemispherical fluxes to the common surface shared by two adjacent cells were summed for the equivalent of a 4π integration.

Discrete ordinates method (DOM)

The discrete ordinates method was first introduced by Chandrasekhar [START_REF] Chandrasekhar | Radiative Transfer[END_REF] and has been widely used since then in radiative transfer applications. A major advantage of this method is that it can be simply coupled with the hydrodynamics system using the same structured or unstructured grid [START_REF] Murthy | [END_REF]. The DOM is based on the discretization of the RTE (see Eq. ( 1)) over a chosen number, N dir , of discrete directions, s i (µ i , η i , ξ i ), contained in the solid angle 4π and associated with weights w i . Koch and Becker [9] compared several types of angular quadrature schemes, of which the two that were selected for this study are S 4 and LC 11 . The S 4 quadrature scheme, with 24 directions, is often chosen for its computational efficiency. The LC 11 scheme, with 96 directions, was chosen for its higher accuracy. The ordinate directions for each of these quadrature schemes are tabulated in Table 1.

The RTE is solved for every discrete direction s i using a finite volume approach. The integration of the RTE over the volume V of an element limited by a surface Σ with outer unit normal n, and application of the divergence theorem yields (6)

Σ I s • n dΣ = V κI b -κI (s) dV .
The domain is discretized in control volumes which, in this study, are regular hexahedra, but can also be a nonregular mesh for this formulation. Taking I j to be the average intensity over the j th face, associated with the center of that face, and taking I b,P and I P to be the average intensities over the volume V , associated with the center of the cell, P , Eq. ( 6) can be discretized as ( 7)

N face j =1 I j (s i • n j )A j = κV (I b,P -I P ),
where the scalar product of the ith discrete direction vector with the normal vector of the j th face of the considered cell is defined by s i • n j = µ i n xj + η i n yj + ξ i n zj . I b is assumed to be constant and equal to I b,P over the volume V , and I j is taken constant over each face. For each cell, the incident radiation G given in Eqs. ( 2) is evaluated at the center by

(8) G ≃ N dir i=1 w i I P (s i ).
For a gray medium, the divergence of the radiative heat flux is obtained from Eqs. (2). To solve Eq. ( 7), a spatial differencing scheme based on the mean flux (DMFS), proposed by Ströhle et al. [START_REF] Ströhle | Int. Conference on Heat Transfer[END_REF], was used. This scheme uses the decomposition (9)

I P = 1 2 I out + 1 2 I in ,
where I in is the weighted average of the intensities at the entering faces of the cell and I out the weighted average of the intensities leaving the cell. Substituting I out from Eq. ( 9) into Eq. ( 7), and after some algebra (see [START_REF] Joseph | [END_REF] for details), the following expression for I P results,

I P = 1 2 κV I b - j,D ij <0 D ij A j I j (10) 1 2 κV + j,D ij >0 D ij A j ,
where D ij = s i • n j . After I P is calculated from Eq. ( 10), the radiation intensities at those cell faces at which D ij > 0 are set equal to I out , obtained from Eq. ( 9). In solving the RTE along a given discrete direction, the control volumes should be treated following a sweeping order, depending on the considered direction and defined such that the radiation intensities are known at the upstream cell faces. A specific algorithm is then used to define the order of the computation of the cells. A radiative heat transfer code called DOMASIUM has been developed based on this method to compute the radiative transfer in complex geometries, mainly for combustion applications [START_REF] Joseph | [END_REF]. Radiation is solved on hybrid grids (hexaedra and tetraedra), imposed by the CFD code developed at CERFACS 3 which solves the hydrodynamics by LES techniques. To take into account gas combustion spectral lines, an SNB-ck model has been implemented. In this 3 Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique.

work, only regular cells and simple gray medium are considered; nevertheless, this comparison represents a validation of the code in its simplest formulation.

Monte Carlo method-net exchange formulation (MCM-NEF)

Monte Carlo methods (MCM) have often been used to produce highly accurate solutions in the process of validating other numerical methods [42,43]. They first appeared in the literature as strict numerical implementations of stochastic photon transport models [44,[START_REF] Hammersley | Monte Carlo Methods[END_REF]. The very large number of realizations required to achieve convergence shows the limitations of the classical Monte Carlo algorithms, particularly when optically thick media are encountered [START_REF] Farmer | Radiative Transfer: Current Research[END_REF][START_REF] Farmer | [END_REF]. To overcome these difficulties, a mathematical formulation using the net exchange formulation (NEF) [48,49], together with adapted probability density functions, has been proposed to improve the variance reduction procedures [50]. In addition, to treat each ray, computer graphics techniques are applied to reduce the computational time. A new space subdivision called voxels, defined below, enables the use of efficient algorithms for intersection calculations [START_REF] Perez | Int. Conference in Central Europe on Computer Graphics, Visualization and Computer Vision[END_REF].

Taking P i as a point within the volume V i and P j within V j , we denote the position vectors of P i and P j as r P i and r P j . The net radiative exchange between two volumes V i and V j , ϕ (V i ,V j ) , or a volume V i and a surface S j , ϕ (V i ,S j ) , or two surfaces S i and S j , ϕ (S i ,S j ) , is expressed for black walls and nonscat-tering media as ( 11)

ϕ (V i ,V j ) = V j V i κ(r P i )κ(r P j )τ (s ij ) s 2 ij × I b (r P i ) -I b (r P j ) dV i dV j , (12)
ϕ (V i ,S j ) = S j V i |n(r P j ) • s|κ(r P i )τ (s ij ) s 2 ij × I b (r P i ) -I b (r P j ) dV i dS j , (13)
ϕ (S i ,S j ) = S j S i |n(r P i ) • s||n(r P j ) • s|τ (s ij ) s 2 ij × I b (r P i ) -I b (r P j ) dS i dS j ,
where s ij = s js i = |r P j -r P i |, and

(14) s = (r P j -r P i ) |r P j -r P i | ,
with n the normal vector to the surface S, κ the gray absorption coefficient, and τ (s ij ) the spectral transmissivity along a straight line between P i and P j given by (15) τ (s ij ) = exp -s j s i κ(s) ds .

The radiative source term for a volume V i and the net heat flux at a surface S i are computed by taking into account their radiative exchanges with all the other volumes and surfaces, (16)

S r (r P i ) = V i ∇ • q r dV i = N s j =1 ϕ (V i ,S j ) + N v j =1 ϕ (V i ,V j ) and (17) q w,net,i = N s j =1 ϕ (S i ,S j ) + N v j =1 ϕ (S i ,V j ) ,
where N s is the number of surfaces and N v the number of volumes.

One way of evaluating the multiple integrals in the expressions for the net exchange rates, Eqs. ( 11), [START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF], and [START_REF] Mihalas | Foundation of Radiation Hydrodynamics[END_REF], is to use a Monte Carlo method, which is described next.

Considering that each radiative exchange can be represented as an integral I of a function f over a domain D,

(18) I = D f (x) dx.
An arbitrary probability density function (PDF), p, defined and strictly positive on the integration domain D is introduced. The weight function W (x) = f (x)/p(x) is used to write

I = D f (x) p(x) p(x) dx = D W (x)p(x) dx.
Given a random variable X, distributed according to p, and a function of that variable, g(X), we let I represent the expectation of g(X). Estimating I with N samples g(x i ), where x i is the ith realization of the random variable X gives

I = E g(X) ≈ 1 N N i=1 g(x i ) = g(X) N , (19) where
I = lim N→∞ g(X) N .
Then the standard deviation of the estimate is calculated as σ (⟨g(X)⟩

N ) = 1 √ N σ (g(X)
), where σ (g(X)) is the standard deviation of g(X). It is approximated by

(20) σ g(X) N ≈ 1 √ N g(X) 2 N -g(X) 2 N .
In this last expression the variance depends on the function g, which itself depends on the PDF. To perform efficient Monte Carlo simulations, the choice of the PDF is crucial. More details can be found in [50,[START_REF] Eymet | [END_REF].

It should be noted that the results presented in this paper have a standard deviation of about 1% or less. Moreover, this Monte Carlo algorithm has another feature that combines the PDF optimization with computer graphics techniques in order to optimize the ray-tracing procedure. It applies an efficient algorithm to evaluate the intersections between paths and surfaces and uses voxels. Voxels, or "volume elements," are cells defined in terms of homogeneous emission temperatures and molar fraction levels that subdivide the space, but differently from the CFD grid [START_REF] Perez | Int. Conference in Central Europe on Computer Graphics, Visualization and Computer Vision[END_REF]. To save computational time, a good compromise between the time to treat each voxel and the level of grid subdivision has to be determined.

Ray tracing

This method treats the RTE as a set of firstorder ordinary differential equations (ODEs), with one ODE for each spatial point, directional angle, and frequency. In the gray problem considered here, the frequency dependence is handled by analytically integrating the RTE over frequency. At each spatial point x, a set of rays is considered to project inward toward the point, with the set being chosen to sample solid angle space in such a way as to allow accurate integration over that space to compute the net incident radiation and the heat flux. For the fire problem, the rays are followed outward from the chosen point until they intercept a wall or the exit of the chimney. At such a boundary point the initial value of the incoming radiative intensity is set to equilibrium (I = I b ). From this initial value the RTE is integrated forward along the ray to the chosen spatial point, and the value at that point is saved for inclusion in angular integrals involving I (x, Ω).

The method of integration along a ray assumes that, within each step of the quadrature, the source I b and opacity κ are constants equal to their interpolated values at the center of the step. With this assumption, I (s, Ω) can be advanced from one end of the step, s 0 , to the other, s 1 , according to the following rule:

I (s 1 , Ω) = I (s 0 , Ω) exp -κ|s 1 -s 0 | (21) + I b 1 -exp -κ|s 1 -s 0 | .
Once the full set of angular values I (x, Ω) is obtained at the point x, angular integrals, such as those in Eqs. (2), are performed to compute quantities of interest.

Moment methods and the M 1 closure

A system of equations for two moments, the net incident radiation G and the radiative flux q r , can be extracted from the gray RTE, Eq. ( 1), by integrating over all directions. The system is given by ( 22)

1 c ∂ t G + ∇ • q r = κ 4σ T 4 -G , (23) 1
c ∂ t q r + ∇ • (D r G) = -κq r ,
using the same notation for the other solution methods above. The M 1 closure [23,24] is given by the Eddington tensor D r obtained by the maximum entropy closure [10,25]. More references concerning this model can be found in [26]. The Eddington tensor is computed from the Eddington factor χ and the anisotropic factor f = q r /G as

(24) D r = 1 -χ 2 Id + 3χ -1 2 f ⊗ f f 2
, where Id denotes the identity matrix, f the Euclidean norm of f, and ⊗ stands for the dyadic product. The Eddington factor χ is given by ( 25)

χ(f ) = 3 + 4f 2 5 + 2 4 -3f 2 .
The Eddington tensor D r , which plays the role of a flux limiter, comes from an underlying radiative intensity that is able to describe both a beam by a Dirac function and isotropic radiation by a Planck function.

Hence, the M 1 model is able to predict radiation in opaque, semiopaque, or transparent media and, as we show below, is particularly suited for the computation of radiation in fires. The numerical scheme used to solve this model is given in [53].

Results and discussion

In the following computations, the Monte Carlo code provides the radiative source term (Eqs. ( 2)) of reference. The ray-tracing code obtains the same values for this term when angular convergence is ensured. Angular convergence of this code is reached for 20,000 angles. This has been checked with comparing with a 80,000 angles computation. The ray-tracing code with 20,000 angles is thus also considered as a second reference code and provides the reference radiative fluxes, which our MC code does not provide (only by a lack of implementation).

Angular resolution and coupling

In this section, the number of angles needed for the baseline fire computation is investigated using the ray-tracing code. This result should aid future choice of angular resolutions, quadrature schemes, and solution methods. Fig. 4 shows the radiative source term, div(q r ), plotted as a function of elevation above the fuel pan for several radial positions with distinct radiant transport behavior. Inside the fire at a radius of r = 0.43 m, high emission temperature radiation is incident from all directions, creating an isotropic environment. At larger radii within the fire (e.g., at r = 0.7 m), conditions are such that the radiation from high-temperature soot is near isotropic. At all positions at high elevations above the fuel pan and at radial positions outside the fire (i.e., r > 1.0 m for this fire), the field is anisotropic since high emission temperatures exist only on one side of the point of interest. The reference solution is shown with 20,000 angles (µ = 100 × φ = 200), which was sufficient for angular convergence. 4 It is apparent that for low angular resolutions, less than approximately 50 angles, predictions of the source term are poor both inside (Figs. 4a and4b) and outside (Figs. 4c and4d) the fire. When only 5 or 10 angles were used, which is likely during the initial stages of a simulation to accelerate it toward a steady solution, it was found that ray effects were dominant, meaning that a hot source may or may not transport heat to each distant point according to the angles chosen. Results vary greatly based on the choice of angles for integration. Inside the fire (Figs. 4a and4b), in the isotropic regions, it is found that at least 50 rays are needed to get results close to the reference profile. Outside the fire (Figs. 4c and4d), it is shown that the use of 350 angles gives good agreement and ray effects are reduced or eliminated. Hence, because such resolutions are needed for accuracy, a high computational cost is expected. Nevertheless, these results must be balanced by the fact that neither special quadratures, nor particular choices of angles have been used herein to try to improve the accuracy of the results for low angular resolution. This feature will be illustrated in the next sections. In the ray-tracing solver, angles are uniformly distributed in µ = cos(θ ) and φ, which may not be the optimal choice. Undoubtedly, a better choice of angles and/or quadratures could decrease the number of angles needed to get accurate results.

Fig. 5 shows profiles of the radial component of the radiative flux as a function of elevation at radial locations inside (Figs. 5a and5b) and outside (Figs. 5c and5d) the fire. Similarly to the source terms, it is apparent that angular resolutions of less than 50 rays lead to poor prediction of the flux both inside and outside the fire. Moreover, in the higher anisotropic regions outside the fire, more than 250, ideally 350, rays are required to reasonably match the reference solution.

The importance of the ray effects are emphasized by considering a coupled fire problem. Fig. 6 shows the predicted hemispherical heat flux incident on a surface at a fixed location in the fire at radius r = 0.55 m and height z = 0.5 m above the fuel pan as a function of time in the 2-m fire. The source term predicted at each angular resolution was coupled to the hydrodynamics and evolved in time from ignition. Radiation was solved by the DT method for various angular resolutions from as low as 24 rays to as much as 350 rays (as suggested by the results in Figs. 4 and5). Of these results, it is seen that 24 angles does not lead to an accurate solution of the radiative flux and contains strong temporal fluctuations. Since the coupling between radiation and hydrodynamics is strong in fires, a poor computation of the radiation can, for instance, lead to extinction or to an over/under-estimation of soot formation. Since these phenomena are highly nonlinear and strongly coupled, chained effects on the resulting fires can occur as in the two following scenarios. The prediction of excess soot (which corresponds to high opacity) can lead to excessive heat loss by radiation, which then leads to a smaller, weaker fire, or even to extinction. On the other hand, an underestimation of soot can lead to low radiation heat losses, implying a hotter, larger fire. Thus, low angular resolution for radiation should be avoided for fire computations to ensure proper coupling and corresponding fire behavior particularly during transients. It is also seen from Fig. 6 that using 80 rays seems to be a good compromise between speed and accuracy for the Vulcan results to be close to convergence, but the fluctuations are still notable. It is reasonable, then, as a result of these observations, to run a fire simulation with Vulcan with 80 angles to allow the fire to develop more rapidly, and later use 350 angles for a converged solution with less fluctuations. For the remainder of this study, all Vulcan results are reported for the higher, converged solution at a resolution of 350 angles obtained after the initial progression with 80 angles to a steady-like state.

Radiative source term

A comparison of the radiative source term (div(q r )) calculated by each method is shown in Fig. 7. Since this quantity couples the hydrodynamics system to the radiant transport solution, accurate computation is mandatory for obtaining the correct fire profile. A positive value of the source term represents a net emission, and vice versa, negative values imply that net absorption is occurring at that point. The ray-tracing code with high angular resolution and the Monte Carlo code were used to determine reference solutions of this quantity. Both of these codes found similar solutions, within the 1% accuracy of the Monte Carlo technique, at all points in the facility. Figs. 7a and7b show the source term profiles for radial positions r = 0.29 m and r = 0.43 m inside the fire as a function of elevation above the pan. Good global agreement with minimal discrepancies is found between all methods. The differences are sufficiently small that they would not have a strong effect on the coupled energy equation. The similarity in accuracy of these results can be explained by the fact that the radiation inside the fire is mainly isotropic, f < 0.2, which would lessen ray effects.

For points just outside the fire at r = 1.15 m in Fig. 7c, more discrepancies are seen among the methods. The lower accuracy is likely due to the less isotropic nature of radiation in this region; however, it is noted that this far outside the fire the small errors would not affect the dynamics of the fire. The DOM S 4 method is less accurate than DOM LC 11 , most likely because its lower angular resolution of 24 angles (versus 96 of LC 11 ) results in ray effects. The M 1 model gives results nearly identical and as accurate as DOM S 4 . These methods slightly underestimate absorption outside the fire, while the DTM overestimates it to a similar degree. Despite the discrepancies, the magnitude of the source term is sufficiently small in this region that less accuracy is acceptable.

In Fig. 8, two-dimensional contour plots of the source term along the central plane of the fire are in comparison of the solution methods. 6 In this figure, no noticeable difference is apparent in the source term computation among the four methods considered here. In particular, the large positive values, indicating emission and representing the radiative cooling, are the same, as shown in detail in Figs. 7a and7b. This finding is encouraging for the low-cost M 1 and S 4 methods, as accurate determination of the source term ensures a solution coupled to the hydrodynamics that will evolve properly. Negative values of the source term, indicative of absorption zones, are also visible in Fig. 8 inside the fire and indicate the preheating cone of the fire. In this region, the chemistry evolves strongly and most of the radicals and pollutants are formed or preformed. Moreover this zone is heated by absorption of the radiation coming from the hot regions of the fire, which is here correctly computed by all methods. To conclude, it has been found that all six methods give similar predictions of the radiative source term. The fact that S 4 computes both flux and source terms relatively well, when linked with the angular studies of Section 4.1, implies that 24 angles should be enough inside the fire provided that the S 4 set of angles is chosen. This result constitutes an improvement by a factor of two of the number of angles needed, compared to the use of uniformly distributed angles. (It should be also noted that the 24 angles of DTM used in the time evolution study in the previous section were not optimally chosen as are the 24 rays in the S 4 technique.)

Radial component of the radiative flux

Fluxes outside the fire are shown in Figs. 9e and9f. The DTM method gives accurate results, in agreement with the ray-tracing solver. The DOM LC 11 also gives accurate results, which leads to the conclusion that around a hundred angles should be enough to compute radiation at 1 or 2 diameters away from the fire provided an accurate quadrature is chosen. This constitutes an improvement by a factor of 3 compared to a uniformly spaced angular set, which apparently needs 350 angles (see Section 4.1). The DOM S 4 does inside the fuel pan and the wall (the interior zones are here visible, as for instance next to the chimney, in the upper right corner). The immersed boundary conditions, as used here to model this complex geometry, were found to be efficient even though an artificial absorption zone was found in the first cells of the pan in direct contact with the fire. There is no artificial flux created below z = 0 and no heat travels through the walls or the pan. All methods considered found flux vectors similar in direction, independent of the method used; the global direction 1 m from the centerline which is slightly tilted north east is found by the four methods considered here. Moreover, the radiation emitted from the hot regions and absorbed in the preheating cone can be seen. Fig. 10 also illustrates a ray effect in both DOM solutions at an elevation just above the pan at z = 0.1 m, evident by the sharp indentation in the contours of constant flux at radii larger than the pan, r 1 m, and stronger fluxes downward below the pan level. The same effect could be seen previously in the profiles in Figs. 9e-9f. This effect is believed to be due to proximity to the fuel pan. Both DOM methods suffer from this ray effect, although LC 11 is significantly closer to the ray-tracing solution than S 4 , and is thus the more accurate solution technique. On the other hand, the flux contours of M 1 do not contain the same sharp indentation nor strong downward fluxes below the pan elevation. This is believed to be due to the M 1 closure which, by construction, allows only one anisotropic direction at one given location. As a result, a global, unique, and, as previously described, slightly tilted northeast direction is observed. This direction is believed to be due to the hot emitting region at r = 1 m. Along this direction, the fluxes taken at a diameter away from the fire begin to be overestimated by M 1 , as seen previously in Figs. 9e-9f. A nozzle profile is then formed due to the fact that this model does not limit sufficiently the fluxes in regions of strong anisotropy.

Conclusions

Six different methods were used to compute the radiative field of a synthetic 2-m, JP-8 pool fire, including Monte Carlo, ray tracing, DOM S 4 , DOM LC 11 , DTM, and the M 1 moment model. The M 1 method was applied for the first time to a combustion problem occurring in a complex three-dimensional geometry. Theoretically, this model has the lowest computational cost of the six, since the directional integration is handled analytically with a model.

An angular resolution study has shown that a minimum of approximately 50 angles inside and 350 angles outside the fire are needed to accurately compute radiation when a uniformly distributed set of angles is chosen. This choice of angles is not optimum, however, and the number of angles needed can be reduced when an optimum quadrature is chosen. It was shown that these numbers can be reduced to 24 inside when the S 4 quadrature is chosen and to 100 outside when the LC 11 quadrature is chosen. Unfortunately, using this many angles still results in a considerable computational cost. It was also shown that if an insufficient angular distribution is used, significant changes to the time-dependent solution may occur. It is thus not possible to accurately compute a time-accurate fire if the aforementioned angular requirements are not satisfied.

The ray-tracing and Monte Carlo methods, which are the most accurate methods when their convergence is ensured, were used to compute reference solutions for the radiative source term and radiative flux fields. Results of both methods were shown to be essentially identical. Both of these methods, which need a large number of angles or a large realization sample, respectively, are too costly to be used for a threedimensional time-dependent fire. It was found that inside the fire in isotropic regions all six methods give similar predictions of the radiative source term, which is needed for coupling with the hydrodynamics. Close to the fire, this term is underestimated by the S 4 and M 1 methods and overestimated by DTM, but the deviations were small enough that they would not affect the dynamic fire behavior. This demonstrates that the M 1 method gives results within the fire of sufficient accuracy to be recommended for the computation of fire dynamics.

More discrepancies were found in the comparison of the radiative flux. Inside the fire, all methods show agreement due to the isotropic conditions which exist inside the domain. However, outside the fire, where anisotropic conditions exist, both the S 4 and M 1 methods inaccurately compute the flux. DOM S 4 exhibited ray effects and the M 1 method overestimated the maximum flux by almost a factor of 2. These two methods are thus not effective if high accuracy is required for the heat flux to an object far outside the fire. On the other hand, DOM LC 11 and DTM gave accurate solutions outside the fire.

Overall, it was shown that M 1 results are most similar to DOM S 4 in the sense that its results are accurate in regions where S 4 results are accurate, and therefore, it is recommended that the M 1 method be considered for future combustion applications, since it is far less expensive computationally. The discrete transfer method was shown to give results similar to the ray-tracing reference solution, which was expected; however, the reader is reminded that it is not as efficient as DOM S 4 and M 1 . Finally, DOM LC 11 gave results very close to the Monte Carlo and raytracing reference codes and is recommended for situations where accuracy is more important than efficiency, especially in anisotropic regions of the computational domain.

Fig. 1 .

 1 Fig. 1. (Left) The FLAME facility for fire experiments and (right) a still image from video taken during a puff of a 2-m JP-8 pool fire test.

Fig. 2 .

 2 Fig. 2. (Left) A two-dimensional contour view of the computational domain modeling the FLAME facility and (right) the emission temperature (( T 4 ) 1/4 ) contours along the centerline planes of the three-dimensional synthetic fire.

Fig. 3 .

 3 Fig. 3. A two-dimensional contour view of the absorption coefficient κ.

  (1) dI (s) ds = κI b -κI (s), where I b = σ T 4 /π is the blackbody intensity at temperature T , σ is the Stefan-Boltzmann constant, and κ is the filtered absorption coefficient (based on a Planck mean coefficient). The gray intensity I represents here a spectral intensity I (ν) which has been integrated on the whole frequency spectrum ν ∈ [0, +∞[. I b accounts for the temporal fluctuations of the fourth power of the temperature due to TRI, and thus I does too.

Fig. 4 .

 4 Fig. 4. Radiative source term, div(q r ), in (W/m 3 ) as a function of elevation above the fuel pan computed by ray tracing with the specified angular resolution µ × φ at radial positions (a) r = 0.43 m, (b) r = 0.70 m, and (c, d) r = 1.03 m.

Fig. 5 .

 5 Fig. 5. Radial component of the radiative flux (W/m 2 ) as a function of elevation above the fuel pan computed by ray tracing with the specified angular resolution µ × φ at radial positions (a) r = 0.43 m, (b) r = 0.70 m, (c) r = 1.03 m, and (d) r = 2.55 m.

Fig. 6 .

 6 Fig. 6. Temporal evolution of the radial component of the hemispherical heat flux to a surface at an elevation z = 0.5 m and radial position r = 0.55 m with a radially outward normal calculated by the discrete transfer method with various angular resolutions.

Fig. 7 .

 7 Fig. 7. Radiative source term (W/m 3 ) as a function of elevation above the fuel pan at radial positions (a) r = 0.29 m, (b) r = 0.43 m, and (c) r = 1.15 m.

Fig. 8 .

 8 Fig. 8. Contour plots of the radiative source term (W/m 3 ) along the central plane of the fire domain computed by, from left to right, ray tracing, DOM LC 11 , M 1 , and DOM S 4 . Radial and axial distances have units of meters.

Fig. 9

 9 Fig.9shows the radial component of the radiative flux in W/m 2 as a function of elevation above the fuel pan for radial locations inside and outside the fire. The ray-tracing code at a high angular resolution was used to generate the reference solution since the computation of the flux field by the Monte Carlo code was not yet available. In Figs.9a and 9bthe fluxes inside the fire are compared. All methods roughly agree and compute consistent trends. At higher elevations, for z > 1.5 m, all methods fully agree since

Fig. 9 .

 9 Fig. 9. Radial component of the radiative flux (W/m 2 ) as a function of elevation above the fuel pan at radial positions (a) r = 0.29 m, (b) r = 0.43 m, (c) r = 0.70 m, (d) r = 1.03 m, (e) r = 1.42 m, and (f) r = 2.00 m.

Fig. 10 .

 10 Fig. 10. Contour plots of the radial component of the radiative flux (W/m 2 ) along the centerline of the fire domain computed by, from left to right, ray tracing, DOM LC 11 , M 1 , and DOM S 4 .

Table 1

 1 Ordinate direction components of S 4 and LC 11

	Quadrature	µ	η	ξ	w
	S 4	0.2958759	0.2958759	0.9082483	0.5235987
		0.2958759	0.9082483	0.2958759	0.5235987
		0.9082483	0.2958759	0.2958759	0.5235987
	LC 11	0.1891433080	0.1891433080	0.9635609052	4π/96
		0.1891433080	0.9635609052	0.1891433080	4π/96
		0.9635609052	0.1891433080	0.1891433080	4π/96
		0.4486223380	0.4486223380	0.7729657144	4π/96
		0.4486223380	0.7729657144	0.4486223380	4π/96
		0.7729657144	0.4486223380	0.4486223380	4π/96
		0.6865960430	0.6865960430	0.2391061427	4π/96
		0.6865960430	0.2391061427	0.6865960430	4π/96
		0.2391061427	0.6865960430	0.6865960430	4π/96
		0.8795381380	0.4758283974	0.	2π/96
		0.4758283974	0.8795381380	0.	2π/96
		0.8795381380	0.	0.4758283974	2π/96
		0.4758283974	0.	0.8795381380	2π/96
		0.	0.8795381380	0.4758283974	2π/96
		0.	0.4758283974	0.8795381380	2π/96

All ray-tracing results presented were calculated with this resolution.

Monte Carlo and DTM results were not available for every node of the domain and are, therefore, not plotted here.

The very small differences between the emission and absorption terms is amplified by the large opacity chosen for the wall (600 m -1 ).

Acknowledgments

Much of this work was performed during the Center for Turbulence Research 2004 Summer program. The authors express their sincerest gratitude to this institution and, in particular, to its director, Professor Parviz Moin, who gave us the opportunity to accomplish this work. The authors also especially thank Dr. Sheldon Tieszen of Sandia National Laboratories for his extensive technical discussions and encouragement. Kirk Jensen was supported by the Advanced Simulation and Computing Program of Sandia, a multiprogram laboratory operated by Sandia Corporation, a Lockheed-Martin Company, for the United States Department of Energy's National Nuclear Safety Administration under Contract DE-AC04-94AL85000.