Two-step method for radiative transfer calculations in a developing pool fire at the initial stage of its suppression by a water spray

A procedure based on two-step method is suggested to simplify time-consuming spectral radiative transfer calculations in open flames containing scattering particles. At the first step of the problem solution, the P1 approximation is used to calculate the divergence of radiative flux, and it is sufficient to determine the flame parameters. The second step of solution is necessary to obtain the radiation field outside the flame, and this can be made independently using the raytracing procedure and the transport source function determined at the first step. Such a splitting of the complete problem results in much simpler algorithm than those used traditionally. It has been proved in previous papers that the combined two-step method is sufficiently accurate in diverse engineering applications. At the same time, the computational time decreases in about two orders of magnitude as compared with direct methods. An axisymmetric pool fire at the initial stage of fire suppression by a water spray is considered as the case problem. It is shown that evaporating small water droplets characterised by a strong scattering of infrared radiation are mainly located in regions near the upper front of the flame and one can observe the scattered radiation. This effect can be used in probe experiments for partial validation of transient Computational Fluid Dynamics (CFD) simulations.


Introduction
Radiative transfer calculations are the most time-consuming in CFD simulation of fires.
Therefore, the use of a simple but sufficiently accurate method for the radiative transfer is critically important to improve the computational costs of CFD codes. A numerical analysis of suppression of fires by water sprays is possible only with the use of spectral models because the optical properties of water droplets cannot be considered on the basis of a gray model [1]. In most cases, there is no need for a very detailed radiation field in CFD calculations. Only the total or integral (over the spectrum) radiative flux divergence is important because this term is needed in the energy equation that couples thermal radiation with other modes of heat transfer. It means that one can consider a differential approach instead of the radiative transfer equation (RTE). Such a choice is really reasonable in the case of combined heat transfer problems including those for fires and combustion systems [1,2]. The P 1 (the first-order approximation of the spherical harmonics method) [1][2][3][4] is the simplest method of this type, and there is a very positive long-time experience of using this approach in diverse multi-dimensional problems [1,[5][6][7][8][9][10][11][12][13][14][15][16][17]. It appears to be sufficiently accurate in the case when one needs only the divergence of radiative flux in the energy equation. Moreover, it was shown in recent paper [18] that the P 1 error is sometimes less than that of the finite-volume method.
It was demonstrated in early papers [19][20][21] that P 1 approximation may give incorrect values of radiative flux near the boundaries of the computational region in the case of a strong decrease of temperature or extinction coefficient of the medium towards the region boundaries. Unfortunately, this is the case for open flames. Therefore, the second step of the computational radiative transfer model is necessary. This can be done using the transport approximation [12] and a ray-tracing procedure for the RTE with the transport source function determined at the first step of solution. In particular, this combined two-step method presented in papers [5,19,20] was successfully employed in recent paper [17] to calculate radiative heat transfer from supersonic flow with suspended particles to a blunt body.
The present study is focused on applicability of the earlier developed combined two-step method to the flames typical of pool fires. The small developing flame just before its suppression by water sprays and also at the initial probe stage of suppression is considered. It should be noted that radiation scattering by evaporating droplets is important even at the beginning of fire suppression and this effect is favourable for the P 1 accuracy because of more smooth angular dependence of the spectral radiation intensity [1].
Following paper [22], we consider infrared radiation outside the absorption band of water to focus on the effect of radiation scattering by water droplets. The simple box model is used for spectral absorption coefficient of molecular gases, but this choice is not critical for the present study because it is focused on a general approach for the radiative transfer calculations and also on infrared scattering by evaporating water droplets outside the absorption bands of water. A numerical solution for the conventional axisymmetric flame (after averaging the numerical results for the 3-D flow field) is considered because this simplification in the flame data is acceptable for the present study.

An estimate of the radiation effect on the flow field calculations
The effect of radiation on the flow field of a small developing flame before the flame suppression is considered. The fire scenario for the case problem is based on a 28 cm diameter methane gas flame, burning in the open quiescent environment with a heat release rate of 53 kW.
The LES (large eddy simulation) CFD simulation was carried out using an in-house version of FireFOAM code [23], compatible with OpenFOAM version 2.2.x. The combustion sub-model used in this code is based on the Eddy Dissipation Concept (EDC) proposed in paper [24] for RANS (Reynolds-Averaged Navier-Stokes) applications and extended to LES in [25]. The radiative transfer equation is solved using the FVM code, which is open accessible. Gas radiation is treated with the box model approach based on the exponential wide band model to calculate the equivalent band absorption coefficient [26]. The flow field in the present study was computed using FireFoam.
The calculated transient flame is three-dimensional, but deviations of the main parameters from the axisymmetric case are not significant. Therefore, we consider similar axisymmetric fields obtained by an averaging in a selected axial cross section of the real developing flame at = 1s. A comparison of two temperature fields shown in Fig. 1 indicates that the effect of thermal radiation on the flame formation is relatively small at early stages of flame development. In particular, the maximum temperature in the flame is underestimated by about 1.5 % in the calculation without taking into account radiation. Nevertheless, thermal radiation should be taken into account because of its considerable effect on the parameters at the upper "head" of the flame (compare Figs. 1a and 1b). In light of this observation, rather than using detailed and time-consuming methods for radiative transfer calculations, a much simpler P 1 approximation seems to be reasonable. Of course, the error of this simple approach should be additionally verified to be sure that it is acceptable for the developing flame calculations. As to the radiative flux from the flame, it can be calculated using a ray-tracing procedure at the second step of solution independently of the CFD calculations. athe calculation without account for thermal radiation, bthe complete calculation.

Two-step method for radiative transfer calculations
The general theoretical model for transfer of thermal radiation in a flame should be based on the radiative transfer equation (RTE) which takes into account emission, absorption, and scattering of radiation in the medium. A special feature of the fire suppression problem is the anisotropic scattering of radiation by relatively cold water droplets. In this case, the RTE can be written as follows [1] (hereafter, we omit the subscript " " in the medium radiative properties for brevity): The physical meaning of Eq. (1) is evident: variation of the spectral radiation intensity, λ , in direction Ω ⃗⃗⃗ takes place due to self-emission of thermal radiation by gases (the last term, b ( ) is the Planck function), extinction by absorption and also by scattering in other directions (the second term), as well as due to scattering from other directions (the integral term). The spectral extinction coefficient is a sum of the absorption coefficient, = g + d , (both gas and water droplets contribute to absorption) and the scattering coefficient, of a conventional continuous medium of droplets. These spectral coefficients and the scattering phase function, Φ, depend on spatial coordinate ⃗.
The integration of RTE (1) over all values of the solid angle gives the following important expression for divergence of the spectral radiative flux [1]: which can be treated as a local balance of the spectral radiation energy. The spectral radiative flux and irradiation function are defined as follows: The first term in the right-hand side of Eq. (2) is the emission, whereas the second term is the radiation absorption by both the gases and water droplets.
A complete and accurate solution to the RTE in scattering media is a very complicated task. One can find many studies in the literature on specific numerical methods developed to obtain more accurate spatial and angular characteristics of the radiation intensity field. Several modifications of the discrete ordinates method (DOM), finite-volume method (FVM), and statistical Monte Carlo (MC) methods are the most popular tools employed by many authors (see [1][2][3][4] for the details). The main difficulty of accurate radiative transfer calculations is related with the angular dependence of radiation intensity. Fortunately, this angular dependence appears to be rather simple in many important applied problems. It enables one to use this property of solution to derive relatively simple but fairly accurate differential approximations.
All the differential approximations for RTE are based on simple assumptions concerning the angular dependence of the spectral radiation intensity λ . These assumptions enable us to deal with a limited number of functions λ i ( ⃗) instead of function λ ( ⃗, Ω ⃗⃗⃗ ) and turn to the system of the ordinary differential equations by the use of integration of RTE. With the use of the spherical harmonics method, the spectral radiation intensity is presented in a series of spherical functions. In the first approximation of this method, P 1 , the following linear dependence is assumed [1]: By multiplying Eq. (1) by Ω ⃗⃗⃗ and integrating it over a solid angle, one can find that the radiative flux is related to the irradiation by equation is the spectral radiation diffusion coefficient, tr = + tr is the transport is the transport scattering coefficient, ̅ is the asymmetry factor of scattering introduced as follows: Note that Eq. (5) is obtained for an arbitrary scattering function but it is the same as that for the transport approximation when this function is replaced by a sum of the isotropic component and the term describing the peak of forward scattering [1,12]: where 0 = Ω ⃗⃗⃗ . Ω ⃗⃗⃗ ′ is the cosine of the angle of scattering. It means that P 1 approximation is insensitive to details of the scattering function and the asymmetry factor of scattering is the only parameter of scattering anisotropy taken into account in this approach. The latter does not lead to considerable errors in radiative heat transfer calculations in the practically important case of multiple scattering [1,12].
Substituting Eq. (5) into the radiation energy balance (2) we obtain the nonhomogeneous modified Helmholtz equation for the spectral irradiation: The Marshak boundary conditions are used to complete the P 1 boundary-value problem [1]. One can also proceed to the variational formulation of the problem when the solution yields the minimum of the following functional [1,5]: where w and w are the temperature and hemispherical emissivity of boundary surfaces. The first integral in Eq. (9) corresponds to the differential Eq. (8), which is true in the computational region, whereas the second integral takes into account the conditions at the region boundary. The order of the differential operator in Eq. (9) is reduced by one compared with Eq. (8). As a result, the local linear approximations of function λ in triangle finite elements are sufficient. The final-element method (FEM) [27,28] can be used to solve the variational problem for functional (8) as it was done in papers [9,17,29]. The FEM is very convenient because this method is flexible and applicable for computational regions of complex shape.
It should be recalled that the box spectral model with constant values of absorption coefficient in a set of spectral bands is used in the current version of the CFD code FireFOAM. According to this spectral model, the Planck function, b ( ) , in Eq. (8) is replaced by the following integral: As a result, the integral value of irradiation, j , for every band of number "j" is obtained instead of the spectral irradiation. A contribution of the spectral band to the divergence of the integral (over the spectrum) radiative flux is given by the following relation: Let consider the computational data for CO 2 absorption band with wavelength boundaries of 1 = 4.1494μm and 2 = 4.5851μm for the model flame with temperature field shown in Fig. 1b. The absorption coefficient in this spectral band calculated with the box gas radiation model is plotted in Fig. 2. The radial optical thickness,  of the hottest region of the flame without water droplets is easily estimated by multiplying the absorption coefficient, g,j , by the geometrical thickness, Δ , of this thin region. The calculations carried out by the authors show that a value of that is much smaller than unity is obtained. It means that the upper estimate of the radiative heat losses based on the optically thin approximation may be applicable in this scenario specific to the initial development stage of the fire.
The non-uniform triangulation of the computational region, which includes the main part of the temperature field of Fig. 1b, is shown in Fig. 3. The left boundary of the region is the surface of the liquid fuel. The smaller areas of finite elements were made in the region of more optically dense medium. The total number of triangular elements is 2x40x40=3200. As usually, the radial coordinates of all the nodes were increased by the first radial interval of the mesh to avoid the formal difficulties in numerical calculations near the axis [1]. The finite elements near the front of the developing flame should be relatively small (as it is shown in Fig. 3) because of an expected increase in extinction in this part of the flame containing numerous absorbing and scattering water droplets. Note that triangular elements at a small distance from the liquid fuel look not good because of a significant difference between the radial and axial dimensions of the triangles.
However, additional calculations showed that this weakness is not critical in the case of the optically thin developing flame. The following dimensionless function is calculated using various approaches: This function is convenient to compare the results of calculations with those obtained in the optically thin limit when there is no need to solve the equations presented above because the following simplest relation is true for the radiative flux divergence [3,4]: The results of calculations are presented in Fig. 4. The FVM calculations are carried out using FireFoam. One can see that the difference between the computational data obtained using the optically thin approximation, P 1 method, and FVM are similar to each other and the first of these approaches gives an upper estimate of radiative heat losses. The calculations based on P 1 and FVM give very close values of the radiative flux divergence, whereas the computational time of the P 1 calculations is about two orders of magnitude less than that of the regular FVM procedure used in the CFD code FireFoam. The results obtained confirm good accuracy of P 1 approximation which can be employed for very fast and reliable radiation calculations for the developing flames and can be recommended for implementation into any general CFD code. It should be recalled that P 1 is even more accurate in the case of a scattering medium [1] and can be used in radiative heat transfer calculations for the flames with water droplets as applied to the flame suppression problem.
where λ is the so-called source function. The spectral radiation intensity can be calculated by direct integration of Eq. (14). This procedure is very simple in the case of an open flame when no radiation illuminates the flame surface and the formal solution along the ray ⃗ is as follows: Eq. (14) can be rewritten as applied to a set of spectral bands considered in the box model. In the combined two-step method of the present paper, the source function is calculated at the first step using of the P 1 solution for irradiation λ . After that, a ray-tracing procedure for the transport RTE solution should be used at the second step to calculate accurately the radiative flux to an arbitrary surface outside the flame, if necessary. It is interesting that there is a potential possibility of comparison of computational predictions with experimental data for the radiation of pool fire flames [30].

Computational model for the probe stage of fire suppression
The effects of water sprays on the flow field parameters of the flame should be taken into account in computational studies for the regular developed regime of fire suppression, and this is an important engineering problem [31][32][33][34]. However, the problem is significantly simplified at the probe stage of fire suppression suggested in [22] because this preliminary stage of fire suppression is characterized by a very small flow rate of water. The observations at the probe stage could show us some details of behavior of evaporating water droplets in the flame.
A short-time probe stage of fire suppression is considered in the present paper. Therefore, the motion, heating, and evaporation of single droplets can be studied without account for any feedback effects. The radiative power absorbed by droplets is neglected because it is much less than the heat supplied to the droplets by convective heat transfer from hot combustion gases in the flame. In more detailed calculations, the effect of radiation can be taken into account as it was done in recent papers on shielding of fire radiation by water mists [35,36] and also in the problem of a space vehicle shielding from solar radiation by a cloud of sublimating particles [37].
The interaction of water sprays with fires has been modeled computationally in many papers during the last two decades [38][39][40][41][42][43][44][45], but there is no need to discuss here the state-of-the-art in this field because of the main focus on the radiative transfer problem. Of course, there is a size distribution of water droplets at the initial cross section of the spray. However, the assumption of independent behavior of droplets makes it possible to consider a set of solutions for droplets of different size in the expected range of the initial droplet size distribution. The trajectories of droplets moving along the flame axis are considered. The radial component of gas velocity in a flame is much less than the axial gas velocity. In addition, the initial direction of the droplet motion is also close to the axial one. This makes the 1-D approach for the dynamic and thermal interaction of water droplets and flame to be acceptable for the computational estimates of the main effects. To simplify the problem, we do not consider also the secondary effects of the vapor coming from evaporating droplets on both the drag force and convective heat flux to the droplet. This is the usual assumption used in computational models of this kind. The initial problem for the droplet motion is as follows [18,41,46]: where the subscript " " refers to water, is the axial component of velocity, is the droplet radius, D is the drag coefficient, Re is the Reynolds number. It is assumed that water droplets are first heated up to the saturation temperature:  (18) where w is the latent heat of evaporation of water. A similar model has been also used in papers [22,35,36]. This simplified evaporation model is also sufficient for qualitative estimates of the present paper because of very fast heating of moving droplets and large molar fraction of water vapour in the flame. It is not necessary at the moment to consider sophisticated evaporation models like that developed in paper [47]. A transfer from the droplet heating to its evaporation is given by the following simple equation: Obviously, the local relative volume fraction of water droplets can be calculated as: Let consider the axial profiles of a gas temperature, g , and velocity, g , on the axis of Equation (21) gives very small values of rel < 0.25ms for water droplets of radius < 100μm. It means that water droplets should be supplied at a small axial distance from the flame and the initial velocity of droplets should be sufficient to reach a hot region inside the flame.
Obviously, the effect of initial conditions decreases along the droplet trajectory, and the equilibrium position of a droplet is characterized by the same absolute values of two opposite forces: the force of gravity and the drag force which supports the droplet. In the case of a relatively small evaporating droplet, this condition can be written using the Stokes law. The resulting value of the gas axial velocity is as follows: Having substituted It is physically clear that the volume fraction of water droplets increases dramatically in the local regions of their collection and total evaporation. Strictly speaking, the assumption of independent motion of single droplets may be incorrect in the vicinity of these regions. However, following [20], this effect is not considered in the present paper.

Optical properties of water droplets
The normalized coefficients of absorption and transport scattering (per unit volume fraction of water) for independently absorbing and scattering droplets with radius can be determined as follows: According to Mie theory [49][50][51], the efficiency factors of absorption, , and transport scattering, The spectral optical constants of water, and , can be found in early papers [53,54]. According to [55,56], we neglect a weak temperature dependence of the refractive index. It means that the most important value of transport scattering coefficient in the range of water semi-transparency can be calculated using the values of optical constants of water at room temperature. The calculated values of a and s tr presented in Fig. 7 should be multiplied by the volume fraction of droplets to obtain the local values of the spectral absorption and transport scattering coefficients: d = v a and tr = v s tr . One can see that radiation absorption by water droplets is significant in the main nearinfrared absorption band of water. It contributes to water evaporation and flame suppression. The scattering by small droplets is large over the spectrum excluding the absorption band, where the anomalous dispersion takes place [57,58]. The hypothesis of independent scattering does not work in the region of strong evaporation where the distance between the neighboring droplets may be comparable with the wavelength. The dependent scattering may lead to significant changes in optical properties of the medium [59][60][61][62]. However, this effect is a complex problem of physical optics and it is beyond the scope of the present paper. Therefore, the ordinary relations of independent scattering by single droplets are used in the estimates. It should be noted that there are two serious assumptions of the suggested theoretical model and both of them are related with an expected high local volume fraction of small water droplets.
Strictly speaking, one cannot consider independent motion of these droplets and also independent scattering of infrared radiation by water droplets in some local regions of the flame. The resulting overall error of the simplified computational model depends strongly on water flow rate at the probe stage of the flame suppression. One can decrease this error in the experiments using additional probe-stage regimes with much lower flow rate of water.

Effect of radiation scattering by evaporating water droplets
The computational data presented in Fig 7 confirm that the above choice of the representative CO 2 absorption band at wavelength about 4.3μm will enable us to estimate the effect of scattering of the flame radiation by evaporating water droplets at the probe stage of fire suppression. Note that contributions of soot to both absorption and scattering of radiation are negligible in this spectral range [63][64][65][66]. The absorption of thermal radiation by droplets is very small in this range whereas the normalize transport scattering coefficient is very high: s tr = 300cm −1 at = 10μm and s tr ≈ 1000cm −1 at = 5μm (see Fig. 7b). Following [67], the emissivity of the liquid fuel surface w = 0.6 was used in the calculations. The latter does not mean that thermal radiation of relatively cold liquid fuel is important. However, we take into account a partial reflection of the flame radiation from the surface of liquid fuel.
It was shown in the previous section of the paper that relatively small evaporating droplets of It is assumed that scattering by water droplets in the region of < cl and also at > 0.2m is negligibly small.
As discussed above, the observed radiation from the flame is determined by the radiative source function. Therefore, the effect of infrared scattering by water droplets on the source function is considered. Following the above normalization, the source function for the transport RTE calculated using the P 1 solution for irradiation is divided by the maximum value of the source function calculated using the optically thin approximation. The finite-element triangulation of the computational region shown in Fig. 5 was used in the calculations. The resulting field of dimensionless value of S calculated using P 1 approximation is presented in Fig. 8. One can see that a cloud of small evaporating droplets leads to significant changes in the field of radiative source function just in front of the developing flame. This effect is expected to be interesting for the infrared diagnostics of transient flames at the probe stage of fire suppression. The latter can be used for partial validation of CFD simulations. To the best of our knowledge, there is no appropriate data in the literature on infrared measurements of thermal radiation of a flame suppressed by a jet of water droplets. However, such an experimental work would be interesting.

Conclusions
A procedure based on the earlier developed combined two-step method is suggested to simplify time-consuming spectral radiative transfer calculations in transient pool fires during their suppression using water sprays. A comparison with the reference numerical solution indicates that the accuracy of the P 1 approximation used at the first step of the method is sufficient for the calculations of the radiative flux divergence. This will lead to a significant simplification of the radiative module in the CFD code and much faster radiation calculations. This result can be considered as an important stage in improving the fire radiation modeling.
The calculations of motion, heating, and evaporation of water droplets at the probe stage of fire suppression indicated the local areas of focusing of moving droplets before their total evaporation in the flame. The effect of infrared radiation scattering by these small droplets on the source function responsible for the radiation from the flame is considerable. It is important that positions of the local highly-scattering areas are very sensitive to the flame parameters. This can be used in probe experiments for partial validation of transient CFD simulations.

Conflict of interests
None declared.