Cyril Caliot
email: cyril.caliot@promes.cnrs.fr

Vincent Eymet

Mouna El-Hafi

Yannick Le Maoult

G Flamant

M El Hafi

Parametric study of radiative heat transfer in participating gas-solid flows

Keywords: Monte Carlo, Non-gray gases, Spherical particles, Anisotropic scattering, Refractive index

come

Introduction

Thermal radiation, in high temperature gas-solid flows, is considered to play a significant role with respect to the other heat transfer mode. Two-phase mixture is common in engineering practice, e.g. particle-gas heat exchanger, pulverized coal furnace, solid rocket motor etc. But the presence of various radiating gases, soot particles and polydispersed particulates complicates significantly the study of such processes. For example, combustion processes generate gaseous species, mostly CO 2 and H 2 O, and particles (soot, fly-ash particles or carbon aggregates).

In such situations, modelling of thermal radiation is far from being complete due to many complexities such as radiative properties of non-gray gases, multi-dimensional geometry and computational accuracy as well as efficiency. In addition particles absorb, emit and scatter anisotropically thermal radiation. Their radiative properties are strongly dependent on composition, which affects their optical properties. Moreover, the particle size distributions and the volume fraction affect the extinction coefficient, the scattering albedo and the angular distribution of the scattered intensity of the particulate phase. For these reasons, modelling radiative transfer in a high temperature gas-solid flow is computer time consuming. To reduce this computation time, the prediction of the radiant exchange uses various methods, models and assumptions depending on the nature of the problem to solve. A better understanding of the major radiative processes involved in the studied application can be reached by performing a sensitivity analysis. As an example, in order to study potential simplifications in engineering applications of coal furnaces, isotropic and neglecting of scattering are evaluated against anisotropic scattering in [START_REF] Marakis | A parametric study of radiative heat transfer in pulverised coal furnaces[END_REF]. Furthermore, isotropic and anisotropic scattering heat transfer processes have been analyzed in a parametric study [START_REF] Yuen | A parametric study of radiative transfer with anisotropic scattering in a one-dimensional system[END_REF]. Another parametric study was performed in [START_REF] Yu | An extension of the weighted sum of gray gases non-gray gas radiation model to a two-phase mixture of non-gray gas with particles[END_REF] to examine the non-gray radiative effects of the two-phase mixture in an axisymmetric cylinder by varying parameters such as the temperature of particles, the non-gray gas composition and the particle concentration.

At the beginning of a survey and before building a radiative model for a gas-solid flow, one must have identified the main radiative transfer contributions involved. Then, respective contributions of gas and particles to the whole emission and the influence of scattering have to be evaluated. Gas spectra have a fine structure that is produced by line transitions. Consequently, the spectrum has to be divided into narrow bands and the treatment of the radiative transfer involves great resources in computer time and memory. For particles, the multi-scattering processes complicate the radiative model and require additional computing resources. From that point of view, it is recommended to construct a radiative transfer model in a gas-solid flow after firstly estimating the role of gas emission/absorption and particles scattering.

This paper aims to present a sensitivity study which permits, on the one hand, to evaluate the influence of scattering on the radiant exchanges in a particular gas-solid flow and, on the other hand, to assess the influence of gas emission/absorption in such a flow. For this purpose, a Monte-Carlo method is used to compute the radiant heat flux leaving a one-dimensional homogeneous layer containing non-gray gases (H 2 O and CO 2 ) and a monodispersion of spherical particles which scatter anisotropically. Since it is difficult to measure the optical properties of particles in a two-phase flow, this paper addresses this difficulty by considering a wide range of refractive index and by performing a sensitivity analysis.

In this paper, the radiative transfer model using the Monte-Carlo method and the CK spectral model is presented first. Secondly, the methods for the computation of the radiative properties of gases and particles are given. Then, a sensitivity study is reported to show the dependency of particle radiative properties on their optical properties. Moreover, the influence of scattering by particles and gases radiation on the spectrum emitted by the layer is emphasized.

Description of the radiative transfer model and the radiative properties

To reach a better understanding of the influence of scattering as well as the influence of gas emission/absorption in a realistic gas-solid flow, an approach is proposed, in this paper, which consists to realize parametric studies on simplified situations.

The first simplification addresses the 3D geometry involved in realistic problems. The complex 3D geometry of a specific application can be simplified to a 3D slab (1D geometry). This simplification is chosen to consider only one geometrical parameter: the layer thickness e (Fig. 1); which is characteristic of the realistic problem. Obviously, the thicknesses of the two other dimensions are infinite. Then, the influence of both scattering and gas contribution on the emitted flux by the layer (Φ in Fig. 1) is maximized since optical thicknesses are maximum in the perpendicular directions of the slab thickness e. The second simplification addresses the distribution, into the slab, of particles and gases concentrations, as well as temperature. In this study, a homogeneous slab with constant concentrations and temperature is considered. This simplification is done to perform reference sensitivity studies of several homogeneous media which are intended to be compared to realistic and non-homogeneous media. Then, the sensitivity results give one some information about the influence of scattering and gas radiation on the emitted flux by a realistic media. But, prior to the use of the sensitivity results, one has to decide how to use them. If the situation is close to a homogeneous one, the sensitivity results deliver quantitative informations. But, if the homogeneous assumption is far from being realistic, one can divide the medium into layers and apply the sensitivity results to find quantitative informations for each homogeneous layer.

In the following, the study deals with analysis of the emitted radiation from a 1D homogeneous layer containing a two-phase mixture (Fig. 1). The layer contains H 2 O and CO 2 gases with a monodispersion of spherical particles at high temperature. The computation of emitted fluxes by the 1D layer is achieved by a Monte-Carlo Method (MCM) allowing one to solve radiative transfer including phenomena of emission, absorption and anisotropic scattering.

The Monte-Carlo simulation code that was used in the present paper is named KARINE. 1 The underlying algorithm is presented in [START_REF] Eymet | A boundary-based net exchange Monte-Carlo method for absorbing and scattering thick media[END_REF]; its convergence qualities have been improved on the basis of methodological developments that consisted in both a reformulation of radiative transfer in terms of Net Exchanges [START_REF] Cherkaoui | Radiative net exchange formulation within 1D gas enclosures with reflective surfaces[END_REF][START_REF] Cherkaoui | Monte-Carlo simulation of radiation in gases with a narrow-band model and a net-exchange formulation[END_REF][START_REF] Dufresne | Méthode de Monte-Carlo par échanges pour le calcul des bilans radiatifs au sein d'une cavité 2D remplie de gaz[END_REF][START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward Monte Carlo methods and a correlatedk approach[END_REF] and an optimization of sampling laws [START_REF] Cherkaoui | Monte-Carlo simulation of radiation in gases with a narrow-band model and a net-exchange formulation[END_REF][START_REF] Martin | Monte Carlo analysis of the backscattering of radiation from a sphere to a plane[END_REF][START_REF] Dufresne | Inverse Gaussian k-distributions[END_REF]. These efforts permitted to solve the well known convergence problem of traditional Monte-Carlo algorithms in the case of strong optical thickness configurations [START_REF] De Lataillade | A net exchange Monte Carlo approach to radiation in optically thick systems[END_REF]. The formulation used in [START_REF] De Lataillade | A net exchange Monte Carlo approach to radiation in optically thick systems[END_REF] has been generalized in order to take into account the scattering phenomena. The algorithm is fully 3D, but it is used in KARINE for simulations of radiative transfer in monodimensional configurations [START_REF] Eymet | Simulation Monte-Carlo et analyse en puissances nettes échangées des transferts radiatifs infrarouge avec diffusion : vers une paramétrisation dans un modèle de circulation générale atmosphérique[END_REF], such as planetary atmospheres.

This algorithm allow one to compute the total flux emitted by the homogeneous slab. Furthermore, the algorithm takes into account the emission and absorption of the two-phase mixture along a 3D multi-scattering path with the length ℓ. The total flux emitted by the two-phase layer can be written as:

Φ = 1 2 ∞ 0 dλ V dV 4π dΩ(⃗ s) Γ (V,S) P Γ (γ ) dγ κ λ × exp[-κ λ ℓ γ ]L bλ (T ) (1)
with dV an elementary volume in the layer and dΩ the elementary solid angle in the ⃗ s direction (dΩ = sin(θ ) dθ dφ with θ and φ respectively the altitudinal and azimuthal angle). κ λ is the total monochromatic absorption coefficient (particles and gases) and L bλ is the monochromatic black body intensity at the temperature T of the layer. One can notice that the definition of dV and κ λ implies there is enough particles in dV to assume the radiative properties to be homogeneous. Γ (V , S) is the whole optical paths emitted by V and reaching the S surface (layer's border). γ is the variable of the optical path including multi-scattering over all Γ possibles (γ ∈ Γ ). P is the probability density function of appearance of γ . Every multi-scattering paths are computed with a random step model. ℓ γ is the length of the optical path with multi-scattering.

The spectral model used for the two-phase mixture is a k-distribution model with the Correlated-K (CK) assumption [START_REF] Caliot | Modélisation et simulation de l'émission énergétique et spectrale d'un jet réactif composé de gaz et de particules à haute température issus de la combustion d'un objet pyrotechnique[END_REF][START_REF] Taine | Gas IR radiative properties: from spectroscopic data to approximate models[END_REF]. The spectral treatment with the MCM algorithm is detailed in [START_REF] Eymet | A boundary-based net exchange Monte-Carlo method for absorbing and scattering thick media[END_REF]. A Gauss-Legendre quadrature of 10 points is chosen for the CK method.

Regarding gases, radiative properties are obtained directly from high resolution spectra computed by a line-by-line method over the spectral range (800-8000 cm -1 ). This spectral range is chosen since we consider that it covers the significant part of gaseous spectra at high temperature. Line-by-line spectra are computed from spectroscopic parameters of CO 2 and H 2 O respectively included in CDSD1000 [START_REF] Tashkun | CDSD-1000, the high-temperature carbon dioxide spectroscopic databank[END_REF] and in the ONERA database [START_REF] Perez | A spectroscopic database for water vapor adapted to spectral properties at high temperature, and moderate resolution[END_REF]. The line by line model and the method employed to compute CK parameters for the gas mixture are presented by the present authors in [START_REF] Caliot | Remote sensing of high temperature H 2 O-CO 2 -CO mixture with a correlated k-distribution fictitious gas method and the single mixture gas assumption[END_REF]. Then, the CK parameters are found for the mixture of H 2 O and CO 2 (CK model with the single mixture gas assumption, [START_REF] Caliot | Remote sensing of high temperature H 2 O-CO 2 -CO mixture with a correlated k-distribution fictitious gas method and the single mixture gas assumption[END_REF]). In addition, the CK method does not produce any error on the emitted flux by the slab due to the correlation problem, since the slab is homogeneous.

Assuming particles as spherical, their radiative properties are computed following the Lorentz-Mie theory [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF], in the spectral range (800-25 000 cm -1 ). Absorption and scattering coefficients of particles, respectively κ p and σ p , are then obtained in addition with the phase function, p, and the asymmetry factor g. The MCM algorithm needs the asymmetry factor of particles to compute the smooth phase function of Henyey-Greenstein, p HG , in order to optimise random samples and save computer time:

p HG (µ) = 1 4π 1 -g 2 [1 + g 2 -2gµ] 3/2
(2)

g = 1 -1 µp(µ) dµ (3)
with µ the cosine of the scattering angle. Radiative properties of particles are considered constant in a narrow band. These radia-tive properties depend on the spectral complex refractive index, m λ = n λik λ , of particles (refractive index of gases is unity, m λ = n λ = 1), the wavelength λ of the incident beam, the diameter D and the volume fraction f v of particles. Considering a monodispersion of spherical particles, absorption and scattering coefficients, respectively κ pλ and σ pλ , are defined by:

κ pλ = 3 2 f v Q aλ (m, D) D (4)
σ pλ = 3 2 f v Q sλ (m, D) D (5)
where, Q aλ and Q sλ are respectively the spectral absorption and scattering efficiency which are computed following the Lorentz-Mie theory. In this paper, absorption and scattering coefficients are computed for a wide range of the parameters such as f v , m (n = ℜ{m}, k = ℑ{m}), D and λ. Radiative properties of particles are computed for a wide range of the particle diameter, 10 -7 < D < 2.10 -4 m. The volume fraction f v is varying between, a weak value encountered in applications, 10 -7 , and a large value, 10 -3 , which is below the limit volume fraction where both dependent absorption and dependent scattering occur significantly (limit f v = 6.10 -3 [START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF]). For convenience, the complex refractive index is considered constant over the spectrum. Two values of the real part of the complex refractive index is chosen, n = 1.7 and n = 2.5. Numerous particles have their refractive index, n, between this two values of ℜ{m}. The particles with n = 1.7 are typically metal oxide and the particles with n = 2.5 are carbon particles. Furthermore, the imaginary part of the refractive index is varying, 10 -4 < k < 20, to address the effects of slight changes in particle composition which affect the optical properties of particles. Two parametric studies are presented in the paper. A sensitivity study on the radiative properties of particles is detailed first. Secondly, the influence of scattering and gas radiation on the radiant flux emitted by various 1D layers is presented.

Sensitivity studies

As described in the previous section, the Lorentz-Mie theory is applied to compute the radiative properties of particles. The corresponding results are analysed in Section 3.1. The Monte-Carlo algorithm is then used to compare the radiant fluxes emitted by scattering and non-scattering 1D layers (Section 3.2). In the same manner, the Monte-Carlo algorithm is used in Section 3.3 to determine the influence of gas radiation. In Section 3.4, a discussion and additional results are provided to extend the findings of the study to other configurations. Finally, the influence of gases and scattering on spectra will be shown in Section 3.5.

Radiative properties of particles

Emission, absorption and scattering of light by a small spherical particle depend on diameter, wavelength and refractive index. In the following, radiative properties of a particle are studied on the basis of a sensitivity analysis that uses diameter and refractive index as parameters. In order to compare several radiative properties of particles at high temperature, the spectral behaviour of their properties is taken into account by considering average quantity over the spectrum. Average values of the absorption and scattering coefficient (from Eqs. ( 4) and ( 5)) are then obtained, in the case of a homogeneous monodispersion, with:

κ * p = ∞ 0 κ pλ L bλ (T ) dλ ∞ 0 L bλ (T ) dλ = 3 2 f v 1 D ∞ 0 Q aλ L bλ (T ) dλ ∞ 0 L bλ (T ) dλ = 3 2 f v Q * a D (6)
σ * p = ∞ 0 σ pλ L bλ (T ) dλ ∞ 0 L bλ (T ) dλ = 3 2 f v 1 D ∞ 0 Q sλ L bλ (T ) dλ ∞ 0 L bλ (T ) dλ = 3 2 f v Q * s D (7)
These mean values, denoted with a star (Eqs. ( 6) and ( 7)), depend on the particle temperature and also on the energy distribution over the spectrum described by the Planck's law. Mappings of Although monochromatic radiative properties are identical for this two temperatures, Planck mean values are greater for T = 2000 K. Indeed, as the temperature increases and according to Wien's law, much importance is given to small wavelength by the Planck's law where the absorption coefficient is large. In Figs. 2.3 and 2.4 the same comment can be done for particles with n = 2.5.

Considering the same particle temperature but with different refractive index, such as in Figs. 2.2 and 2.4, we find the maximum values of Q * a D of particles with n = 1.7 to be higher than those of the particles having n = 2.5. Maxima are reached for optical index values near unity (0.1 < k < 10) [START_REF] Kattawar | Electromagnetic scattering from absorbing spheres[END_REF] and for size parameter x = πD λ close to unity. As the Lorentz-Mie theory predicts [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF], interactions between particles and light are maximum for x = 1. At 2000 K, significant fluxes are emitted in the range (0.5-10 µm) and the diameter range, matching x = 1, is about (0.16-3 µm). This is the same interval of diameter we find in Figs.

Q * s D
values are also maximum for x close to unity. The maximum values of Q * s D are higher for n = 2.5 than for n = 1.7. Furthermore, in the diameter range (0.5-5 µm), where the scattering properties are maximum, a decrease in magnitude occurs for k close to unity (0.5 < k < 1). For the temperature dependence similar comments than those of the previous section can be given: the maximum values of the averaged radiative properties are higher at 2000 K than at 1000 K.

Influence of scattering

The influence of scattering is investigated in a parametric study on the radiative fluxes leaving a 1D homogeneous layer containing gases and particles, both at the same high temperature, 2000 K (Fig. 1). The gaseous phase contains non-gray H 2 O, CO 2 gases and a non-participating gas, N 2 , at atmospheric pressure (10 5 Pa). Molar fraction of the mixture are chosen representative of a fossil fuel combustion, such as x H 2 O = 0.2, x CO 2 = 0.1 and x N 2 = 0.7. The monodispersion of spherical particles is composed of particles with D = 1 µm.

The influence of scattering is studied using a ratio R s which indicates the error due to neglecting scattering in the radiative transfer model:

R s = Φ ns -Φ Φ 100 (8)
where, Φ ns is the total flux leaving the layer which is computed neglecting the scattering processes and Φ is the total flux taking into account scattering (Eq. ( 1)). Nevertheless, for values of k near unity, when absorption is great, R s values are minimum since absorption hide scattering processes. Indeed, as scattering increases the path length of photon bundles, their absorption along this multi-scattering path is increased as well. So, if absorption is large, non-scattered photon bundles or multi-scattered ones, are both attenuated and the emitted flux is no longer affected by scattering processes.

Following the comparison between Figs 4.1 and 4.2, R s ratios are greater for the layer containing particles with refractive index n = 2.5 than the ones with n = 1.7. Moreover, scattering influence appears to be smaller for particles with n = 1.7 than for particles with n = 2.5, specially when optical index k is lower than 10 -2 . This is due to the influence of the scattering coefficient. Indeed, as mentioned in the study of

Q * s D
, scattering for particles with n = 2.5 is greater than scattering for particles with n = 1.7. Figs. 5.1 and 5.2 represent mappings of the R s ratio in the case of a 1 m thickness layer (e = 1 m). Comments on these Figs. 5.1 and 5.2 are similar to the previous one since only the layer thickness is increased to 1 meter. Then, the scattering optical thickness of the layer is greater than the previous one (e = 0.1 m) and the influence of scattering is enhanced. So, neglecting scattering is acceptable in layers with lower particle volume fractions.

In this section, sensitivity studies have shown, for the different test cases described before, that the influence of scattering can be neglected (R s < 5%) for any values of refractive indexes if f v < 6.10 -6 in a 0.1 m thickness layer or f v < 6.10 -7 in a 1 m thickness layer.

Influence of gases radiation

The importance of gases radiation is represented by the use of mappings of the R g ratio:

R g = Φ -Φ p Φ 100 (9)
where, Φ p is the total flux emitted by the layer containing only particles. The standard deviation of the statistical results Φ and Φ p , obtained by a Monte-Carlo simulation, are both lower than 3%. Fig. 6 shows maps of R g versus the volume fraction of particles (D = 1 µm) and the optical index (k) for two values of the thickness layer. The real part of the refractive index is n = 1.7, and the temperature of the mixture is T = 2000 K. Maps in Fig. 6 indicate that gas radiation may be neglected only if the volume fraction is high enough and for k values around unity (0.1 < k < 10). One can notice the shape of the maps which is linked to the shape of 2.2, D = 1 µm). Indeed, when the absorption optical thickness is dominated by particle absorption, gas radiation is absorbed and the emission is principally due to particles. Concerning the results of R g ratio for layers with particle refractive index n = 2.5, they are not represented since results are identical as those obtained for n = 1.7.

Q * a D versus k (Fig.
This section treats the influence of gas radiation, and one can define the limit value of the R g ratio equal to 5% in order to neglect gas radiation. This limit is reached when the medium has a volume fraction f v > 5.10 -6 (e = 0.1 m) or f v > 10 -6 (e = 1 m), both for k values close to unity. In these regions of the kf v domain gas radiation can then be neglected since particle absorption is dominant. For lower k values both gas and particle play a significant role in the emission of the mixture.

Extension of sensitivity results

The extension of the sensitivity results to other situations can not be done quantitatively with the mappings presented in this paper. For such a purpose, other sensitivity studies involving specific parameters have to be considered. But, the present study should be helpful for qualitative studies of slabs involving temperature between 1000 K and 2000 K and particles with refractive index between 1.7 and 2.5. For example, let us consider two slabs at a temperature of 2000 K containing both gases and particles and having a thickness of e = 1 m. They also contains the same concentrations of gases such as the ones given in the previous studies (x H 2 O = 0.2, x CO 2 = 0.1 and x N 2 = 0.7). The first slab, named Slab1, contains a monodispersion of particles with D = 10 µm and n = 1.7 and the second, named Slab2, contains particles with D = 0.7 µm and n = 2.5. The aim of this example is to show how to qualitatively use the sensitivity results given previously. Indeed, the prediction of the limit volume fraction should be done on the basis of a comparison between the particles radiative properties.

Concerning Slab1 (D = 10 µm and n = 1.7), the radiative properties (Figs. 2.2 and 3.2) of its particles are weaker than those of particles with D = 1 µm (n = 1.7). Furthermore, sen-sitivity results for the slab with D = 1 µm and n = 1.7, yield to the following limit value of volume fraction: if f v < 6 × 10 -7 (Fig. 5.1) scattering is negligible, for all values of k. Then, as a consequence of the qualitative approach, one can predict that the limit volume fraction required to neglect scattering, in Slab1, must be higher than 6 × 10 -7 . Indeed, since the radiative properties of particles of Slab1 are weaker than those of particles with D = 1 µm (n = 1.7), the absorbed and scattered parts of the emitted flux in the Slab1 must be lower than the ones in the reference slab, D = 1 µm (n = 1.7) (Fig. 5.1). In order to give a comprehensive proof, the mapping of R s is computed and shown in Fig. 7.1. The qualitative approach is then verified, since the limit volume fraction for scattering is 10 -5 .

Regarding Slab2 (D = 0.7 µm and n = 2.5), the radiative properties of its particles are greater than those with D = 1 µm (Figs. 2.4 and 3.4). Moreover, the sensitivity results of the reference slab (D = 1 µm and n = 2.5) show that scattering is negligible if f v < 6 × 10 -7 (Fig. 5.2), for all k values. Then, using the qualitative approach, one can predict that the limit volume fraction of Slab2 must be lower than 6 × 10 -7 . As in the previous example, the mapping of R s is computed for Slab2 and the sensitivity results are shown in Fig. 7.2. These results

(1) e = 0.1 m, R s ≈ 0%, R g ≈ 34%.

(2) e = 1 m, R s ≈ 47%, R g ≈ 21%. verify the qualitative approach since the limit value of the volume fraction, in the case of Slab2, is 4 × 10 -7 .

The spectral dependency

Previous sections present discussions related to parametric analysis using average radiative properties over the spectrum.

In this section we propose to evaluate the influence of gas radiation and scattering over the spectrum. For that purpose, spectral fluxes emitted by a 1D layer will be compared in Fig. 8. Then, Fig. 8 represents spectral fluxes leaving two slabs with different thicknesses (e = 0.1 m and e = 1 m) accounting for or neglecting radiation from gases and scattering from particles. Each layer (similarly to those presented previously) contains a gas mixture (x H 2 O = 0.2, x CO 2 = 0.1, x N 2 = 0.7) and a monodispersion of spherical particles (D = 10 -6 , m = 1.7i × 10 -2 , f v = 10 -5 ) both at thermal equilibrium (T = 2000 K) and atmospherical pressure (10 5 Pa). Fig. 8.1 shows the spectral fluxes, Φ, Φ ns , Φ p , Φ g emitted by a one-dimensional slab where the thickness is fixed, e = 0.1 (Φ g is the flux emitted by a slab containing only the gas mixture). One can notice that scattering is negligible. Besides, the R s ratio is zero. But gases radiation is important, since R g ≈ 34%. We can observe the common spectral behaviour of the radiative species such as H 2 O and CO 2 . The absorption and emission bands of gases are located in the spectrum around 1.8 µm (H 2 O band), 2.5 (H 2 O band), 2.7 (H 2 O-CO 2 bands), 4.3 (CO 2 band) and beyond (H 2 O band). Fig. 8.2 represents spectral fluxes for a slab when the thickness is increased, e = 1 m. In this case, the influence of scattering is more important since R s is about 47%. The major discrepancies due to scattering neglect arise between the wavelengths 0.5 and 4 µm where scattering coefficient is important. Also, gases radiation has a significant contribution to total emission with a R g ratio reaching 21%. But, particles scatter and, in the same time, absorb a large part of the radiation emitted by the gas mixture, as it can be seen when Φ g and Φ are compared in Fig. 8.2. Then, the amount of radiation emitted by the gas mixture and the particles cloud inside the slabs are decreased by the scattering effects.

Conclusion and future works

The mean scattering and absorption radiative properties of spherical particles at 1000 K and 2000 K were computed for various particles diameter and imaginary part, k, of the refractive index, whose real part was chosen as n = 1.7 and n = 2.5. For the absorption, largest values are obtained in the diameter range (0.1-20 µm) with 10 -2 < k < 2. Moreover, the largest scattering radiative properties are found, for all k, in the same diameter range (0.1-20 µm). The influence of scattering and gas radiation were also investigated. Then, parametric studies of global fluxes emitted by a 1D slab containing H 2 O, CO 2 and particles (D = 1 µm) at 2000 K was performed versus the volume fraction and the refractive index. In this particular situation, for all the studied values of the refractive index, the influence of scattering can be neglected (R s < 5%) if f v < 6 × 10 -6 in a 0.1 m thickness layer or f v < 6 × 10 -7 in a 1 m thickness layer. Furthermore, gas radiation can be neglected obviously when absorption by particles is dominant, i.e., when optical index k is around unity (0.1 < k < 10). In addition, a qualitative approach is proposed to extend the sensitivity results, provided in the paper, to other situations. Finally, the spectral influence of gases and scattering by particles is discussed.

Since this study yields to quantitative results for gases and particles both at thermal equilibrium, the influence of temperature variations can lead to noticeable change in the radiative heat transfer. Another future works of this study address extension to others particle diameters and the development of a validation code which takes into account the scattering processes in a realistic 3D gas-solid flow.

Fig. 1 .

 1 Fig. 1. Schematic representation of the total radiative flux Φ, leaving the 1D layer composed of H 2 O, CO 2 gases and a homogeneous monodispersion of spherical particles, both at thermal equilibrium.

( 1 )

 1 n = 1.7, T = 1000 K. (2) n = 1.7, T = 2000 K. (3) n = 2.5, T = 1000 K. (4) n = 2.5, T = 2000 K.

Fig. 2 .

 2 Fig. 2. Mappings of Q * a D 10 -4 [m -1 ] versus the particle diameter and the optical index k (m = nik), for two values of (n = 1.7, 2.5), and for two values of the particle temperature, 1000 K and 2000 K.

D

  (unity m -1 ) versus the particle diameter and the optical index k (ℑ{m}, m = nik) are shown respectively in Figs.2 and 3.Maps ofQ * a D ,in Figs. 2.1 and 2.2, are obtained for n = 1.7 and for temperature T = 1000 K and T = 2000 K respectively.

  2.2 and 2.4 where radiative properties are maximum. Concerning Q * s D maps, they are shown in Fig. 3 and

( 1 )

 1 n = 1.7, T = 1000 K. (2) n = 1.7, T = 2000 K. (3) n = 2.5, T = 1000 K. (4) n = 2.5, T = 2000 K.

Fig. 3 .

 3 Fig. 3. Mappings of Q * s D 10 -4 [m -1 ] versus the particle diameter and the optical index k, for two refractive indexes (n = 1.7, 2.5), and for two values of the particle temperature.

Fig. 4

 4 represents mappings of the R s ratio for particles with different refractive index, n = 1.7 in Fig. 4.1, or n = 2.5 in Fig. 4.2. The layer thickness is e = 0.1 m (Fig. 1). The standard deviation of the statistical results Φ and Φ ns , computed with the MCM algorithm, are both lower than 3%. Figs. 4.1 and 4.2 show maximum values of R s for high levels of volume fractions. As a matter of fact, influence of scattering is significant when the scattering optical thickness (t * s = σ * p e) is large. One can notice the behaviour of R s with the optical index, k, that is close to the behaviour of Q * s D (Figs. 3.2 and 3.4, D = 1 µm).

( 1 )Fig. 4 .Fig. 5 .

 145 Fig. 4. Mappings of R s [%] versus the volume fraction and the optical index k, for two values of the refractive index, n = 1.7, 2.5 (e = 0.1 m, T = 2000 K).

( 1 )

 1 e = 0.1 m, D = 1 µm, n = 1.7.

( 2 )

 2 e = 1 m, D = 1 µm, n = 1.7.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Mappings of R g [%] versus the volume fraction and the optical index k, for two values of the layer's thickness, e = 0.1 m and e = 1 m, containing particles with refractive index, n = 1.7 (T = 2000 K).

Fig. 8 .

 8 Fig. 8. Influence of scattering and gases radiation on spectral fluxes emitted by two layers having different thicknesses and same composition.

http://www.energetique.ups-tlse.fr/~eymet/karine.html.


David Joseph
email: joseph@enstimac.fr

Mouna El Hafi
email: elhafi@enstimac.fr

Richard Fournier

Bénédicte Cuenot
email: cuenot@cerfacs.fr

Comparison of three spatial differencing schemes in discrete ordinates method using three-dimensional unstructured meshes

Keywords: Radiative transfer, Infrared line spectra, Discrete ordinates, Unstructured grids, Three-dimensional

Comparison of three spatial differencing schemes in discrete ordinates method using three-dimensional unstructured meshes D Joseph, Mouna El-Hafi, Richard A Fournier, B

Introduction

In computational fluid dynamics (CFD), the coupling between radiative heat transfer and combustion is based on the resolution of the energy equation. The heat source term due to radiation is evaluated by taking into account the temperature and radiating species concentration profiles, which are obtained from the solution of the aerothermochemistry equations. Among all the numerical methods developed to calculate the radiative heat transfer, the finite volume method (FVM) and the discrete ordinates method (DOM) offer good compromises between accuracy and computational requirements. These approaches have been widely used to solve radiative transfer problems in structured three-dimensional geometries using Cartesian or cylindrical coordinates. In particular, the DOM, described by Chandrasekhar in 1950 [START_REF] Chandrasekhar | Radiative Transfer[END_REF], has been deeply studied by Lathrop and Carlson in 60-70's [START_REF] Carlson | Transport theory-The method of discrete ordinates[END_REF] and by Truelove, Fiveland and Jamaluddin in the 80's [START_REF] Fiveland | Three-dimensional radiative heat transfer solutions by the discrete ordinates method[END_REF][START_REF] Fiveland | Three-dimensional spectral radiative heat transfer solutions by the discrete ordinates method[END_REF][START_REF] Fiveland | Comparison of discrete ordinates method formulations for radiative heat transfer in multidimensional geometries[END_REF][START_REF] Truelove | Three-dimensional radiation in absorbing-emittingscattering in using the discrete ordinates approximation[END_REF][START_REF] Jamaluddin | Discrete ordinates solution of radiation transfer equation in nonaxisymmetric cylindrical enclosures[END_REF]. Significant improvements have been achieved in the last decade aiming at the reduction of the ray effects and false scattering, more accurate quadratures and the extension to complex geometries. Nevertheless, the coupling between radiative transfer and other physical phenomena, such as combustion and fluid flow at high temperatures, requires the solution of the radiative transfer equation using the same grid employed to solve the other governing equations. Unstructured grids are often used in CFD owing to their geometrical flexibility. In this way, a lot of work has been developed during the last decade to apply the DOM-FVM to non-orthogonal structured grids and unstructured grids in three-dimensional enclosures [START_REF] Raithby | A finite volume method for predicting a radiant heat transfer in enclosures with participating media[END_REF][START_REF] Chai | Finite volume radiative heat transfer procedure for irregular geometries[END_REF][START_REF] Murthy | Finite volume method for radiative heat transfer using unstructured meshes[END_REF][START_REF] Moder | An unstructured-grid radiative heat transfer module for the national combustion code[END_REF][START_REF] Liu | Development of an unstructured radiation model applicable for two dimensional planar, axisymmetric and 3-dimensional geometries[END_REF]. In particular, Sakami and co-workers proposed an accurate method for the spatial discretization by taking into account the exponential extinction [START_REF] Sakami | Application of the Discrete Ordinates Method to combined conductive and radiative heat transfer in a two-dimensional complex geometry[END_REF][START_REF] Sakami | Radiative heat transfer in 3dimensional enclosures of complex geometry by using the discrete ordinates method[END_REF][START_REF] Sakami | A new differencing scheme for the Discrete Ordinates Method in complex geometries[END_REF], but it necessitates to perform a heavy preprocessing procedure. Much less sophisticated schemes are also commonly used. Liu et al. [START_REF] Liu | Development of an unstructured radiation model applicable for two dimensional planar, axisymmetric and 3-dimensional geometries[END_REF] have used the "step" scheme, equivalent of the "upwind" scheme in CFD and Ströhle et al. [START_REF] Ströhle | A mean flux discrete ordinates interpolation scheme for general coordinates[END_REF] have proposed the mean flux interpolation scheme. Another type of spatial discretization was also recently introduced, namely the Discrete Ordinates Interpolation Method (DOIM), which does not rely upon an integration of the radiative transfer equation over the control volumes, but rather on the integration of the equation along a line of sight. This method has been introduced by Seo et al. in 1998 [18] and extended to unstructured grids by Cha et al. in 2000 [19]. Lastly, Koo et al. [START_REF] Koo | Comparison of three discrete ordinates methods applied to twodimensional curved geometries[END_REF] have compared three methods applied to two-dimensional curved geometries: the DOIM, Sakami's approach and the discrete ordinates method in orthogonal curvilinear coordinates [START_REF] Vaillon | Radiative heat transfer in orthogonal curvilinear coordinates using the discrete ordinates method[END_REF].

In our study, a computer code using unstructured meshes and based on the modeling of radiative transfer using the DOM has been developed aiming at a future coupling with a combustion finite volume code available at the CERFACS 1in France. The code has been specially written for unstructured grids using tetrahedrical cells, by trying to avoid complex adaptations that are highly time consuming. The constraints associated to combustion applications where the following:

(1) Accurate coupling with a finite volume CFD code requires that the output is the integrated radiative source within each mesh; (2) The resulting computation time must remain acceptable for realistic industrial geometries in combustion; (3) The line spectra of combustion gases must be accurately represented across the whole infrared range.

Then, to treat general combustion situations, even if the main part of this paper is devoted to gray media, all these constraints lead to the following choices:

(1) Gaseous line spectra properties are represented with a Statistical Narrow Band correlated-k model and parallelization is used to simultaneously compute the radiative contribution of each narrow band;

(2) For each correlated-k computation the radiative transfer equation is discretized with a finite volume approach (in order to avoid additional computational efforts when coupling with the finite volume combustion code) and three spatial differencing schemes are considered: "exponential" [START_REF] Sakami | Radiative heat transfer in 3dimensional enclosures of complex geometry by using the discrete ordinates method[END_REF], "step" [START_REF] Liu | Development of an unstructured radiation model applicable for two dimensional planar, axisymmetric and 3-dimensional geometries[END_REF] and "diamond mean flux" [START_REF] Ströhle | A mean flux discrete ordinates interpolation scheme for general coordinates[END_REF].

The performances of the different spatial discretization schemes have already been widely investigated in cases of simple geometries mapped using Cartesian coordinates, but there is a lack of information regarding the performance of these schemes with unstructured grids. In the present study, after some theoretical derivations (Section 2), we analyze the behavior of the three retained schemes on academic configurations (Section 3). Results are compared with accurate solutions that are analytical or produced by the Ray Tracing Method. In this analysis convergence difficulties are identified in the limit of strong absorption optical thicknesses. These difficulties are further analyzed using a parametric study to cover the wide range of absorption optical thicknesses typical of spectral line gaseous radiation. A particular attention is given to the standard finite volume approximation and its convergence difficulties at high absorption optical thicknesses, independently of the spatial differencing scheme. We then consider configurations with real gaseous line spectra for the whole infrared range in order to explore the effect of such convergence difficulties for combustion applications. This last case, with a complete H 2 O line spectrum, allows to draw first conclusions and to comment on the validity range of three spatial differencing schemes for combustion applications.

The radiative transfer equation (RTE)

Mathematical formulation

Considering an absorbing-emitting and non-scattering gray medium, the variation of the radiative intensity along a line of sight can be written as:

dI (s) ds = κI b -κI (s) (1)
where I (s) is the radiative intensity, I b the radiative intensity of the blackbody, and κ the absorption coefficient. Boundary conditions for diffuse surfaces are taken from the relation giving the intensity leaving the wall I w as a function of the blackbody intensity of the wall I b,w and of the incident radiative intensity:

I w (s) = ε w I b,w + 1 -ε w π n.s ′ <0 I w (s ′ )|n.s ′ | dΩ ′ (2)
where ε w is the wall emissivity, n the unit vector normal to the wall and s ′ the direction of propagation of the incident radiation confined within a solid angle dΩ ′ .

Angular discretization

The DOM is based on the discretization of the radiative transfer equation (RTE) according to a chosen number N dir of discrete directions, s i (µ i , η i , ξ i ) associated with their respective weights w i , contained in the solid angle 4π . In this way, different angular discretizations can be used. A recent study carried out by Koch and Becker [START_REF] Koch | Evaluation of the quadrature schemes for the discrete ordinates method[END_REF] compares the efficiency of several types of angular quadratures widely outside the most common ones. In this study, we choose the S N quadrature, which is one of the most popular and, in order to estimate the role of ray effect, some tests are performed with the LC 11 quadrature as recommended by Koch and Becker.

Spatial discretization for unstructured grids

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 Σ , and the application of the divergence theorem yield:

Σ I (s i ).s i .n dΣ = V κI b -κI (s i ) dV (3)
We assume that I b and I (s i ) are constant over the volume V and equal to their respective averaged value: I b,P and I P (s i ). We also assume that the intensities at the faces are constant over each face and I j denotes the intensity at the j th face. For unstructured grids, the domain is discretized in tetrahedra and every other hybrid grid can be under-meshed by tetrahedra. Then, for tetrahedra, Eq. ( 3) is discretized as:

N face =4 j =1
I j (s i ).(s i .n j )A j = κV I b,P -I P (s i )

where n j is the exit unit vector normal to the surface j . The scalar product of the ith discrete direction vector with the normal vector of the j th face of the considered tetrahedron is defined by D ij :

D ij = s i .n j = µ i n xj + η i n yj + ξ i n zj (5)
The discretization of the boundary condition (Eq. ( 2)) is straightforward:

I w = ε w I b,w + 1 -ε w π n.s i <0 w i I (s i )|n.s i | (6)
For each cell, the incident radiation G is evaluated as follows:

G = 4π I (s i ) dΩ ≃ N dir i=1 w i I P (s i ) (7)
For a gray medium, the divergence of the radiative heat flux S r is given by

S r = ∇.Q r = κ(4πI b -G) (8)
and the radiative net heat flux at the wall is calculated by

Q w ≃ ε πI b,w - n.s i <0 w i I (s i )|n.s i | (9)
For the three different spatial differencing schemes used, calculations of I (s i ) are detailed in the next sections.

The exponential scheme

Sakami et al. proposed in 1998 [START_REF] Sakami | Radiative heat transfer in 3dimensional enclosures of complex geometry by using the discrete ordinates method[END_REF][START_REF] Sakami | A new differencing scheme for the Discrete Ordinates Method in complex geometries[END_REF] a sophisticated spatial scheme for three-dimensional cases that takes into account the exponential extinction along the optical pathway into the cell. Considering an optical path t = s Bs A , where s A and s B are respectively the positions at the entry and at the exit of the cell, the solution of the RTE (Eq. ( 1)) can be written as

I (s, s B ) = I (s, s A )e -κt + κ s B s A I b e -κ(s B -s) ds ( 10
)
The integration over the cell pathlength leads to

I (s, s B ) = I (s, s A )e -κt + 1 -e -κt I b (11)
For a direction s i and a final position s B at the exit face, the position s A at the entry and the associated pathlength t are fixed. The intensity I l going through an exit face l is defined as the averaged intensity over the lth face:

I l = 1 A l A l I (s i , s A )e -κt + I b 1 -e -κt dA l (12)
It has been shown in [START_REF] Sakami | Radiative heat transfer in 3dimensional enclosures of complex geometry by using the discrete ordinates method[END_REF] that, knowing τ i = κl max,i , where l max,i stands for the maximum thickness of the cell for the direction s i , for the three different types of events presented in Fig. 1, one can calculate intensity I l at an exit face as follows:

I l = k D ik <0 A lk A l I k χ i + I b (1 -χ i ) (13)
where A l is the exit surface area, A lk stands for the part of surface area of the downstream face l obtained by projection of the upstream face k on the face l according to s i and χ i is a typical coefficient that accounts for the exponential extinction through the cell:

χ i = 1 A l A l e -κt dA l = 2 τ i 1 - 1 -e -τ i τ i (14)
The difficulty lies in the calculation of the geometric maximum thickness l max,i for all the cases presented in the Fig. 1. l max,i can be determined if we know the coordinates of the points A and B. By referring to Eq. ( 4), I P is directly deduced:

I P = I b,P - 1 κV N face =4 j =1 D ij A j I j (15)
Whereas this scheme avoids negative values for all I j at cell faces, it does not guaranty that I P is positive at the center of each cell, in particular when κ becomes very small.

The step scheme

Much less sophisticated schemes are also available that are worth a detailed attention in our specific context. Considering our constraints in terms of computational times, the use of the exponential scheme will indeed be only justified if it insures significantly higher accuracy levels for combustion applications. In 2000 [START_REF] Liu | Development of an unstructured radiation model applicable for two dimensional planar, axisymmetric and 3-dimensional geometries[END_REF], J. Liu used the "step" scheme, which corresponds to the "Upwind" scheme that is commonly used in CFD. This scheme had already been proposed by Chai et al., in 1995 [9], in order to solve the RTE by the FVM for irregular geometries using curvilinear coordinates. In previous studies, this scheme has often been applied to Cartesian structured grids, in order to avoid negative values that could occur with schemes such as the "diamond" one [START_REF] Chai | Evaluation of spatial differencing practices for the discrete ordinates method[END_REF]. Omitting the scattering phenomenon, the intensity I P is evaluated at the center of the cell by applying:

I P = κV I b - N face j =1 D ij <0 D ij A j I j κV + N face j =1 D ij >0 D ij A j (16)
where N face is the number of faces for the cell (equal to 4 in our case). The downstream surface intensities are set equal to the intensity I P .

The mean flux scheme

In a similar way, Ströhle et al. [START_REF] Ströhle | A mean flux discrete ordinates interpolation scheme for general coordinates[END_REF] proposed a simple spatial differencing scheme based on the mean flux scheme that can be very useful in the case of unstructured grids. This scheme relies on the following formulation:

I P = αI out + (1 -α)I in ( 17
)
where I in and I out are the cell face averaged intensities entering (white cell faces in Fig. 2) and leaving (gray cell face in Fig. 2) the control volume, respectively:

I in = j D ij <0 D ij A j I j j D ij <0 D ij A j (18)
and

I out = j D ij >0 D ij A j I j j D ij >0 D ij A j (19)
Substituting I in from Eq. ( 17) into Eq. ( 4) yields, after some algebra:

I P = αV κI b -Θ j D ij <0 D ij A j I j ακV + j D ij >0 D ij A j ( 20
)
with

Θ = α -(1 -α) j D ij >0 D ij A j j D ij <0 D ij A j (21)
As shown in Fig. 2, if we consider that the scalar product D ij is the transformation term of the projection of the surface A j , following the vector s i , on the plane ∆ orthogonal to this last vector, we can write:

j D ij >0 D ij A j = - j D ij <0 D ij A j = A ∆ (22)
Then for all α in ]0, 1], we have Θ = 1 and finally, Eq. ( 20) takes the following form:

I P = αV κI b - j D ij <0 D ij A j I j ακV + j D ij >0 D ij A j (23)
The case α = 1 corresponds to the Step scheme (Eq. ( 16)) used by Liu et al. [START_REF] Liu | Development of an unstructured radiation model applicable for two dimensional planar, axisymmetric and 3-dimensional geometries[END_REF]. The case α = 0.5 will be called diamond mean flux scheme (DMFS). It can be compared to the diamond scheme used in structured grids but it had been shown that the DMFS used in structured grid is different from the standard diamond scheme because of the mean flux approximation [START_REF] Ströhle | A mean flux discrete ordinates interpolation scheme for general coordinates[END_REF]. It is formally more accurate than the step scheme. After calculation of I P from Eq. ( 20), the radiation intensities at cell faces such that D ij > 0 are set equal to I out , obtained from Eq. [START_REF] Ströhle | A mean flux discrete ordinates interpolation scheme for general coordinates[END_REF]. It can occur that this last term calculated is negative. I out negative is equivalent to the following conditions:

I b,P < (1 -α) - A ∆ κV .I in (24)

Sweeping order optimization

For a given discrete direction, each plane face of each cell is either an upstream face or a downstream face (a face parallel to the considered discrete direction plays no role). The control volumes should be treated following a sweeping order such as the radiation intensities at upstream cell faces are known. This order depends on the discrete direction under consideration. An algorithm for the optimization of the sweeping order has been implemented in the present work. To order the cells, we have to define first the direction s i . To proceed to the ordering, the grid is swept a first time in an arbitrary order with checking, for each cell, if the number of upstream faces where the intensity is known is the same as the number of upstream faces. If this condition is satisfied for one cell, the number of this cell is stored in a list for the direction s i . Then, the intensity of the downstream cell faces can be known, as the upstream cell faces of the next cells which are updated. We should notice that this sweeping order stored for each discrete direction, only depends on the chosen grid and the angular quadrature. That means that the use of different physical parameters do not change the established list of the sweeping order. A sweeping order optimization avoids too many iterations not only for the cases with scattering media and/or reflective walls geometries, but also avoids iteration in the case of black walled enclosures without scattering.

Results and discussion

To compare the efficiency of the spatial differencing schemes for different types of enclosures containing a participating and homogeneous medium, several test cases have been carried out.

It should be pointed out that for output purposes an interpolation is performed on the results obtained with the unstructured code, in order to have the values for a given axis or a given point. This 3D-interpolation can generate inaccuracies, especially for the points taken near the walls or when the grid is coarse.

Application to a black walled cylindrical enclosure

This test concerns the cylindrical geometry to illustrate the influence of the angular quadrature on the radiative heat flux at the walls.

A cylinder (h = 3 m and R = 0.5 m) containing a gray isothermal medium at T = 1200 K is considered. The walls are black and at T w = 300 K. The radiative heat source S r along the axis of the cylinder and the radiative net heat flux Q w at the side wall are obtained with the unstructured code using the three different spatial discretization schemes described previously and the S 8 quadrature is employed. Two unstructured grids are used in this test case (see Fig. 3). A coarse one comprises 4000 tetrahedra and a finer one 210 000 tetrahedra. A third comparison is made on the same fine grid using a more accurate angular discretization (S 12 ). Results are compared to those obtained with the ray tracing method using 320 000 rays. In the case of a homogeneous isothermal medium, the results of the ray tracing method are independent of the grid.

Two different values for the absorption coefficient of the medium have been chosen to represent weak (κ = 0.1 m -1 ) and strong (κ = 10.0 m -1 ) optical thicknesses. It should be noticed that the optical thickness of the cells also depends on the grid refinement. In the case of optically thin media (see Fig. 4), the peak of the net heat flux at the wall predicted by the DOM is about 10% lower than the ray tracing solu- tion. This is due to the fact that the distribution of optical thicknesses in the medium is function of the angular direction and consequently needs a good angular representation (what can be linked to the well known ray effect). Increasing the order of the angular quadrature (Fig. 4(c)) improves the accuracy of the results. However, such an order of accuracy is compatible with today's expectations in combustion problems. For stronger optical thicknesses (see Fig. 5), the ray effect is softened. The aforesaid distribution of optical thicknesses is less sensitive to angular variation. For optically thick media, the increase of the absorption coefficient of the medium yields also an increase of the wall heat flux. The walls receive the radiation mainly from the closest cells of the medium, so there is no need to refine the grid and the angular discretization. In the case of optically thin media (see Fig. 4), the radiative heat source term solutions are in good agreement. Refining the grid gives better results. Increasing the absorption coefficient of the medium leads to decrease the accuracy obtained for the radiative source term solutions (see Fig. 5(a)). The relative difference between the DOM and the ray tracing solutions is negligible in the vicinity of the walls and becomes significant far from the walls. Most of the energy emitted by the medium is absorbed within a short distance. That leads to a strong exponential extinction that none of the three spatial differencing schemes succeeded in representing if the grids are not sufficiently fine (Figs. 4(a) and 5(a)). So, in the next section, a particular attention will be paid to the accuracy on the radiative source term estimation in a simpler configuration. The step scheme gives better results in optically thin case and the DMFS is the most accurate for very absorbing medium. The angular refinement does not really govern the accuracy of the source term along the centerline. We have noticed that the grid refinement influence is stronger than the angular one when we compare the three spatial differencing schemes.

Maximal, mean and minimal relative errors for the source term given in Tables 1 and2 are computed for the three different spatial differencing schemes as follows E = S r,DOM -S r,RT S r,RT × 100 [START_REF] Liu | Application of the statistical narrow-band correlated-k method to non-grey gas radiation in CO 2 -H 2 O mixtures: Approximate treatments of overlapping bands[END_REF] and for the wall heat flux:

E = Q w,DOM -Q w,RT Q w,RT × 100 ( 26
)
with RT standing for the ray tracing. E Mean is simply the averaged values of relative errors calculated at several points. For an intermediate value of κ (κ = 1.0 m -1 ), the errors are also provided in the two tables.

For the wall heat flux, the solutions computed with the three numerical schemes are in a good agreement. Nevertheless, we can notice that the "Diamond Mean Flux Scheme" gives significantly better results than the exponential scheme for κ = 10.0 m -1 . The DMFS remains accurate enough for a large range of optical thicknesses.

Accuracy levels for combustion applications: The high absorption limit under the finite volume approximation

The preceding academic test cases have pointed out two difficulties commonly encountered with such numerical approaches: the ray effect and the convergence difficulties at the optically thick absorption limit. Ray effects are very much dependent on the geometry of the system and can only be reduced by the increase of the number of discrete directions. As far as the optically thick absorption limit is concerned, two remarks can be drawn that motivate further discussions in the present section:

(1) One would expect that the sophisticated exponential scheme lead to better convergence quality than the simple step scheme and the DMFS; which is not the case in the black walled cylindrical enclosures case. of the general meaning of DOM-FVM computations for combustion applications.

In the present section, we argue that the convergence difficulties at the high absorption limit are intrinsically associated with the finite volume approximation and that, for reasons related to spectral correlation effects, accurate radiative heat source fields can be produced with such methods for combustion applications despite of extreme absorption encountered in CO 2 and H 2 O lines centers, in meter scale configurations.

Black walled spherical enclosure: A parametric study of the optical thickness influence

To analyze in more details the trends identified here above, this new section presents a simple test case which provides an analytical solution (of the source term only) in order to further understand, via a parametric study, the influence of the optical thickness on the source term calculations. In this test case, a sphere with a radius R = 1 m has been considered and two different grids have been used (Fig. 6): a coarse one (18 920 cells) and very fine one (140 010 cells).

In a first case, the medium is isothermal at T max = 1200 K and the wall is black at T min = 300 K. In a second test case, wall conditions are identical but the temperature of the medium is depending on the space variable r as follows:

T (r) = T 4 max 1 -r 2 + T 4 min r 2 1/4 ( 27
)
The source term S r,C is evaluated at the center of the sphere using a S 12 quadrature with the three different spatial differ- The relative error of the radiative source term is represented versus κ in Fig. 7 and is computed as follows:

E = S r,C,DOM -S r,C,analytical S r,C,analytical × 100 ( 29
)
From Fig. 7, we can notice that the relative error increases when the absorption coefficient increases. The change of direction of the curve, in Fig. 7(b) and (d), is only due to the use of the absolute value in the relative error formulation (see Eq. ( 29)). In the isothermal test case, extremely poor levels of accuracy are observed when increasing optical thickness (Fig. 7(a) and (c)). That phenomenon could be interpreted as the "false scattering" effect (numerical diffusion linked to cell optical thickness) with the use of the finite volume approximation. This observation can be related to the fact that, in numerical experiments, the geometrical mesh structure is kept constant when increasing the absorption coefficient, leading to a regular increase of the optical thickness of each mesh. For the non-isothermal case, the same phenomenon Absorption coefficient

κ = 0.1 m -1 κ = 1.0 m -1 κ = 10.0 m -1
Relative errors

E Max E Mean E Min E Max E Mean E Min E Max E Mean E Min
Coarse mesh and S 8

Step scheme 0.56% 0.096% 0.0024%

Absorption coefficient κ = 0.1 m -1 κ = 1.0 m -1 κ = 10.0 m -1
Relative errors

E Max E Mean E Min E Max E Mean E Min E Max E Mean E Min
Coarse mesh and S 8

Step scheme 11.0% occurs but radiative exchanges at small distances (exchanges between closest gas volumes) contribute more than the exchanges with the walls. This reduces the optical thicknesses effectively linked to the most contributing exchanges. Consequently, false scattering has less influence and the relative error on the radiative source term is smaller than the one obtained in the isothermal test case. These errors remain however very high when optical thickness reaches values of the order of ten or above, which is commonly encountered at the center of gaseous absorption lines such as those of H 2 O and CO 2 in combustion.

Finite volume approaches in the optically thick absorption limit

The preceding parametric study has indicated that strong errors are to be expected when increasing the cell absorption optical thickness, whatever the numerical scheme among the three considered ones. This seems in contradiction with the fact that the exponential scheme solves the RTE exactly within each tetrahedrical cell in 3D, or triangular cell in 2D. However, when analyzing this point in detail, it appears that the error is entirely associated with the standard finite vol- ume approximation that the intensity is uniform across an upstream face.

This may be seen in the following simple bidimensional case (Fig. 8): an isothermal gray medium of thickness h and absorption coefficient κ (the blackbody intensity being I b ), between two parallel isothermal black walls of surface S (the blackbody intensity being I b,w ). Let us assume that we want to estimate the integrated intensity I f,total impinging at the surface in a direction ⃗ s i orthogonal to the two parallel surfaces. The analytical solution is:

I f,total = S I f dS = I f .S (30)
with

I f = I b,w e -κh + κ h 0 I b e -κr dr = I b,w e -τ + I b 1 -e -τ (31)
where τ = κh is the absorption optical thickness. Whatever the considered numerical scheme, this geometry requires that the field is divided in at least two triangular cells (C1) and (C2), as indicated in Fig. 8. The numerical resolution will therefore include the computation of the averaged intensity I int at the interface between (C1) and (C2). The analytical expression of I int can be easily shown to be:

I int = I b,w χ + I b (1 -χ) (32)
where χ for a two-dimensional system becomes [START_REF] Sakami | Application of the Discrete Ordinates Method to combined conductive and radiative heat transfer in a two-dimensional complex geometry[END_REF]:

χ = 1 -e -τ /τ (33)
Note that those expressions corresponds to Eqs. ( 13) and ( 14) which indicates that the exponential scheme computes I int exactly from the knowledge of a uniform intensity I b,w at the upstream face of (C1), and a uniform emission within the cell. And indeed, the upstream intensity is uniform (emission by the isothermal black surface) and the volume emission is uniform (isothermal medium). At this stage, everything remains therefore exact with a sophisticated enough numerical scheme. However, things are very different if we now consider (C2) to compute the averaged downstream intensity I f from the knowledge of the averaged upstream intensity I int . As no further information is available, the common finite volume approximation consists in assuming that the intensity is distributed uniformly along the upstream face. Under this assumption, the analytical expression of I f may be shown to be (again, this is what the exponential scheme gives for a uniform intensity at the upstream face as it is exact within each cell):

I f = I int χ + I b (1 -χ) = I b,w χ 2 + I b 1 -χ 2 (34)
This expression is approximate. It only matches Eq. (31) at the optically thin limit. Therefore whatever the quality of the spatial integration scheme, an error will occur at the optically thick limit that is due to the lack of information concerning the distribution of the intensity along the upstream face. For illustration, at the two points A and B reported on the Fig. 8 What these figures tell us when looking at the curves labeled "exponential" is the accuracy loss due to the finite volume approximation itself. In the present configuration this is exactly what the exponential scheme would give because the exponential scheme is exact within each cell. We have also reported in the same figures solution and relative errors (as in Eq. ( 29)) resulting of the finite volume approximation combined with the simple approximate step and diamond mean flux schemes. The first comment is that the convergence difficulties at high optical thicknesses are mainly due to the finite volume approximation itself and not to the quality of the numerical schemes. Concerning, the emission and self-absorption modelling, the finite volume approximation is justified when the medium becomes very thick and the error on I f decreases (Fig. 9(b) and (d)). However at such optical thicknesses, the DMFS gives unrealistic negative values for I int which leads to an increasing error on the final intensity and none of the three numerical schemes succeed in modelling the transmission (very small values for I f ).

As a conclusion, we can notice from Fig. 9(a) and (c) that:

• as long as the DMFS does not calculate negative intensities (for I int ), it gives the most accurate approximate solution, • for intermediate optical thicknesses (10 -1 -10 +1 ), the finite volume approximation introduces important inaccuracies in the different numerical schemes (see Fig. 9(a) and (c)).

Accuracy levels associated to finite volume approximations for line spectra of combustion gases

We have seen that even with high accuracy or exact spatial integration schemes, the usual finite volume approximation itself introduces high uncertainties at the limit of optically thick absorption. A last test case is considered here, which is close to combustion applications, in order to show how the problems encountered at high optical thicknesses (at the centers of gaseous spectral lines) may affect the accuracy of the radiative source term in a real combustion gas. The same unity sphere as above is meshed with 30 000 tetrahedra and filled with an isothermal gas containing water vapor and nitrogen. The spectral dependency of the absorption coefficient is here represented with the use of a SNB-ck model [START_REF] Liu | Band lumping strategy for radiation heat transfer calculations using a Narrowband model[END_REF][START_REF] Liu | Application of the statistical narrow-band correlated-k method to low-resolution spectral intensity and radiative heat transfer calculations-Effects of the quadrature scheme[END_REF][START_REF] Liu | Application of the statistical narrow-band correlated-k method to non-grey gas radiation in CO 2 -H 2 O mixtures: Approximate treatments of overlapping bands[END_REF]. The SNB data have been provided by the EM2C laboratory [START_REF] Soufiani | High temperature gas radiative property parameters of statistical narrow band model for H 2 O, CO 2 and CO and correlated-K model for H 2 O and CO 2[END_REF]. The source term S r is obtained by computing N bands × N quad gray calculations where N bands = 367 is the number of narrow bands which have the same spectral width, and for each narrow band N quad = 5 is the number of Gauss-Legendre quadrature points used in the approximation of the spectral integration. Then, the analytical solution integrated over the whole spectrum is written as:

S r,C,analytic = N band i=1 N quad j =1 4π∆ν i w ij κ ij e -κ ij R × (I b,C,ij -I b,w,ij ) (35)
Consequently, using the DOM, the source term [START_REF] Coelho | Benchmark numerical solutions for radiative heat transfer in two-dimensional nongray sooting media[END_REF] is obtained as follows:

S r,DOM = N band i=1 N quad j =1 ∆ν i w ij κ ij (4π I b,ij -G ij ) (36)
where G ij is obtained from Eq. ( 7) for κ = κ ij .

To study the effect of the angular quadrature on the estimation of S r,C , the relative error has been computed according to Eq. ( 29) in a first configuration with X H 2 O = 0.2 and X N 2 = 0.8. S 4 , S 8 , S 12 and LC 11 [START_REF] Koch | Evaluation of the quadrature schemes for the discrete ordinates method[END_REF] have been tested using the DMFS scheme and for each quadrature, the relative error is about 1.2%.

Consequently, the S 4 quadrature is chosen and the relative error of the source term is computed for the three numerical schemes (with the parameter X H 2 O varying). Fig. 10 represents the relative error versus the molar concentration of water vapor. As we noticed previously, the relevant parameter that influences the solution accuracy is the optical thickness which is directly related to the water vapor molar fraction X H 2 O . The sensitivity of the solution accuracy to the molar fraction of water vapor is illustrated in Fig. 10(a). For small molar fractions (X = 0.01), the exponential scheme could not provide an acceptable physical solution because of the sign of the intensity I p which becomes negative. A very small discrepancy between the three schemes is observed, with a relatively better accuracy of the DMFS. Altogether, the accuracy level is much better than what could be expected on the basis of preceding analysis. As the medium is isothermal, the radiative exchanges with the walls are the only relevant ones. They occur at a long distance, which means that the exchanges corresponding to low absorption coefficients contribute more than the one at high values where extinction rapidly takes place. Finally, the last figure Fig. 10(b) deals with a very thick medium where soot particles are added to water vapor and nitrogen, taking into account a very high soot volumetric fraction f v = 10 -5 . The total absorption coefficient is expressed by with [START_REF] Lee | Optical constants of soot in hydrocarbon flames[END_REF]:

κ ij = κ ij,gas + κ i,soot (37)
κ i,soot = 5.5f v ν i (38)
where f v stands for the soot volumetric fraction and ν i is the wave number taken at the center of a narrow band. As soot radiation becomes dominant, the medium tends to become thick and gray. Then, the results obtained in Fig. 10(b) show big discrepancies between the three schemes which confirm the trends of the relative error for S r observed in the first cylinder test case and a good behavior of the DMFS. In this case, because of the gray behavior of soot in each narrow band, each absorption coefficient is increased at the same level and leads to optically thick media. All the frequencies contribute similarly to the estimation of the exchange with the wall, consequently we can use similar explanations than those presented in the gray case. Altogether, these results only confirm the good behavior of the code giving satisfactory accuracy levels for well mixed combustion chambers, where the gas volume can be considered as quasi-isothermal, with intermediate soot concentration levels. It has been shown in the previous section that, for non-isothermal media, the main exchanges contributing to the radiative source term evaluation are the exchanges with the gas at short or intermediate distance. Therefore, a final interesting test case considering the spherical enclosure containing non-isothermal and non-gray medium is necessary but is not presented here because no analytical solution along the radius is available. A Monte Carlo calculation is presently being performed, that will serve as a reference solution for an extension of the present analysis.

Concluding remarks

A DOM code using unstructured grids has been developed for coupling with finite volume CFD-combustion code and has been validated using academic test cases. This study emphasizes the influence of the absorption optical thickness on the accuracy level. We could then show that the difficulties encountered at high optical thicknesses were associated with the standard finite volume approximation where the intensity is distributed uniformly along each upstream face. These difficulties are therefore quite independent of the sophistication level of the spatial differencing scheme. In particular, for such applications as combustion ones, where the constraints in term of computational times are strong, it appears that there is no particular needs for a very accurate but complex scheme such as the exponential scheme, considering that a simple scheme such as the DMFS insures comparable accuracy levels up to intermediate optical thicknesses and encounter the same strong difficulties at the high absorption limit.

These strong difficulties could lead to believe that standard finite volume approximations are useless in the combustion context where strong absorption commonly appears at the center of infrared CO 2 and H 2 O absorption lines. Parametric test cases were considered to check this point and it appeared that when real gas spectra are considered for integrated radiation computations, much less difficulties are met, in terms of accuracy, than what could be expected on the basis of the conclusions of the gray medium analysis.

This conclusion could only be drawn for well mixed combustion chambers and the cases with strongly heterogeneous temperature fields require further analysis.

For such heterogeneous cases, significant short distance radiative exchanges may appear in optically thick spectral regions. As soon as such exchanges play a significant role the standard finite volume approximation will lead to strong errors. It will therefore be required that some informations are kept concerning the distribution of the intensity along the upstream faces. These ideas have been explored by several authors with the finite element approach where the intensity is stored at each nodes (at the cell vertices), which is the better start basis for the description of its distribution along the face [START_REF] Adams | Subcell balance methods for radiative transfer on arbitrary grids[END_REF][START_REF] Wareing | Discontinuous finite element S N methods on three-dimensional unstructured grids[END_REF][START_REF] Dedner | An adaptative higher order method for solving the radiation transport equation on unstructured grids[END_REF]. These techniques have been shown to be very efficient at the limit of strong optical thicknesses. It would be worth looking at their behavior in the limit of high absorp-tion optical thicknesses for combustion applications. However, in the short term, considering today's computational constraints in this field, such approaches could only be considered for computational accurate reference solutions but not for coupling with finite volume CFD-combustion codes. In the long term, it would certainly be of great interest to try and combine Sakami's exponential scheme with such finite element approaches.

NomenclatureA

  surface area . . . . . . . . . . . . . . . . . . . . . . . . . . . . m 2 A ∆ surface area of cell orthogonally planned following s i . . . . . . . . . . . . . . . . . . . . . . . . . . . . m 2 D ij scalar product of s i by n j G incident radiation . . . . . . . . . . . . . . . . . . . W•m -2 I radiation intensity . . . . . . . . . . . . . . . . W•m -2 •sr N dir number of discrete directions N face number of faces of cells Q w net heat flux at the wall . . . . . . . . . . . . . . W•m -2 R radius . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m S r radiative source term . . . . . . . . . . . . . . . . W•m -3 T temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K V volume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m 3 I in averaged intensity over entries . . . . . W•m -2 •sr I out averaged intensity over exit faces . . W•m -2 •sr Q r radiative heat flux vector . . . . . . . . . . . . W•m -2 n unit vector normal to the face s discrete direction unit vector f v soot volumetric fraction h height . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m l max maximum thickness of a cell . . . . . . . . . . . . . . m s coordinate along direction s t optical pathlength through a cell . . . . . . . . . . . m w weight associated to a discrete direction Greek symbols α weighting factor for mean flux scheme ε emissivity κ absorption coefficient . . . . . . . . . . . . . . . . . . m -1 µ, η, ξ director cosines of the discrete direction ν wave number . . . . . . . . . . . . . . . . . . . . . . . . . m -1 Σ surface area delimiting a volume . . . . . . . . . m 2 τ absorption optical thickness Ω solid angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . sr Subscripts A be ginningpointofapathwaythroughacell B endingpointofapathw aythroughacell b blackbody i associated to ith discrete direction of the angular quadrature j associated to j th face of the cell k associated to kth entry of the cell l associated to lth exit face of the cell P associated to the volume of the cell w wall x, y, z Cartesian coordinates

Fig. 1 .

 1 Fig. 1. Six ways to cross a tetrahedrical cell from A to B according to the number of entries and exit faces.

Fig. 2 .

 2 Fig. 2. Geometrical transformation linked to the term D ij .

Fig. 3 .

 3 Fig. 3. Cylindrical enclosure grid.

Fig. 4 .

 4 Fig. 4. Q w on the side wall and S r on the central axis for κ = 0.1 m -1 of a cylindrical enclosure with participating gray medium: (a) Unstructured coarse grid (18 920 cells) and angular quadrature S 8 ; (b) Unstructured fine grid (140 010 cells) and angular quadrature S 8 ; (c) Unstructured fine grid (140 010 cells) and angular quadrature S 12 .

( 2 )Fig. 5 .

 25 Fig. 5. Q w on the side wall and S r on the central axis for κ = 10.0: (a) Unstructured coarse grid (18 920 cells) and angular quadrature S 8 ; (b) Unstructured fine grid (140 010 cells) and angular quadrature S 8 ; (c) Unstructured fine grid (140 010 cells) and angular quadrature S 12 .

Fig. 6 .

 6 Fig. 6. Spherical enclosure grid.

Fig. 7 .

 7 Fig.7. Relative Error for S r,C at the center of the sphere: (a) For the isothermal gray medium case using the coarse grid; (b) For the non-isothermal gray medium case using the coarse grid; (c) For the isothermal gray medium case using the fine grid; (d) For the non-isothermal gray medium case using the fine grid.

Fig. 8 .

 8 Fig. 8. Simplified configuration.

  , the intensity is arbitrarily set to I int = I b,w χ + I b (1χ) whereas its exact value is respectively I A = I b,w and I B = I b,w e -τ + I b (1e -τ ). The resulting solutions and associated errors are represented in Fig. 9 as function of the absorption optical thickness τ for two test cases: Case 1 (transmission): I b = 0 et I b,w = 1; Case 2 (emission/self-absorption): I b = 1 et I b,w = 0.

Fig. 9 .

 9 Fig. 9. Comparison between the four expressions of I f (k) and the relative error for the three differencing schemes: (a) Transmission modelling: I f (k); (b) Transmission modelling: Relative error; (c) Emission/self-absorption modelling: I f (k); (d) Emission/self-absorption modelling: Relative error.

Fig. 10 .

 10 Fig. 10. Relative Error for S r,C at the center of the sphere using the medium grid: (a) Case of the isothermal gaseous medium; (b) Case of the isothermal mixture of gas and soot.

Table 1 Relative

 1 Errors on the radiative source term S r

Table 2 Relative

 2 Errors on the radiative heat flux at the wall Q w

	3.78%	2.60%	0.51%	328%	219%	7.67%

Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique.

Acknowledgements

The financial support of the European Commission in the framework of the TMR programme (RADIARE project, No. ERBFMRX-CT98-224) is gratefully acknowledged. Special thanks for the Pr. Coelho of the IST in Lisboa for his indications and fruitful discussions.