Fatmir Asllanaj
email: fatmir.asllanaj@univ-lorraine.fr

Sylvain Contassot-Vivier

Olivier Botella

Francis H R França

Numerical solutions of radiative heat transfer in combustion systems using a parallel modified discrete ordinates method and several recent formulations of WSGG model

Keywords: Radiative heat transfer, WSGG model, discrete ordinates method, complex 3D geometries, parallel computing, combustion systems. Nomenclature

The aim of the paper is to perform 3D radiative transfer calculations in combustion gas mixtures. For this purpose, a parallel modified Discrete Ordinates Method with a cell-vertex formulation on unstructured tetrahedral grid is proposed for solving the Radiative Transfer Equation (RTE) in realistic geometries. Several recent formulations of the weighted-sum-of-gray-gases (WSGG) model are implemented. To reduce the false scattering, the cell-vertex formulation is combined with an Exponential scheme and is compared to the Step scheme. The results are then compared with published results on non-isothermal and non-homogeneous gas mixtures. The use of non-gray WSGG models shows a very good level of accuracy for the evaluation of wall heat transfer and radiative heat source. The Exponential scheme is more accurate than the Step scheme. The summation over the wavenumber spectrum of the radiative model is usually carried out after having computed the radiative quantities for each discrete angular direction. By summing first over the wavenumber spectrum, a reduction in the computation time is significantly obtained, up to a factor of 16 on a case study. Also, additionnal significant speedup is obtained with an efficient distributed and multi-threaded parallel algorithm. When the Step scheme is used, the computation time increases very slightly with the number of gray gases. The procedure to speed up the calculations can be applied to any arbitrary numerical method for solving the spectral RTE.

Radiative heat transfer plays an important role in many engineering problems. It is often the dominant mode of heat transfer in combustion systems due to the high temperatures (peak temperatures are typically around 2000 K), the presence of radiatively participating species (mainly CO2, H2O, soot, and other particulates in solid fuel flames), and relatively large length scales of such problems [START_REF] Liu | The impact of radiative heat transfer in combustion processes and its modeling -with a focus on turbulent flames[END_REF][START_REF] Coelho | Radiative Transfer in Combustion Systems[END_REF][START_REF] Modest | Radiative Heat Transfer in Turbulent Combustion Systems: Theory and Applications[END_REF][START_REF] Modest | Radiative Heat Transfer[END_REF]. Then, it affects significantly temperature distributions in combustion systems.

Consequently, an accurate evaluation of radiative heat fluxes is crucial if one attempts to improve the efficiencies of such systems, to predict the temperature distribution or to simulate pollutant formation and emission.

Accurate modeling of radiative heat transfer in combustion systems is challenging due to the strong spectral dependence of the radiative properties of combustion products and the huge calculation times required for solving the spectral Radiative Transfer Equation (RTE) [START_REF] Liu | The impact of radiative heat transfer in combustion processes and its modeling -with a focus on turbulent flames[END_REF]. Accurate and efficient RTE solvers are then needed to account for the fine spectral behavior of gas radiative properties. A radiative model is required not only to be acceptably accurate, but also must be computationally efficient to be used in coupled numerical simulation of turbulent combustion problems. This is because radiative heat transfer is only one of the many phenomena to be modeled in turbulent combustion, and the CPU time spent on radiative heat transfer must be kept as small as possible to make the overall modeling task tractable. This is why we have designed in this work a parallel algorithm to solve the RTE.

Most studies on parallel extensions to RTE solvers consider three possible decompositions according to the respective discretizations over angular, spatial and spectral domains [START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF][START_REF] Novo | Parallelization of the Discrete Transfer Method[END_REF]. Concerning the angular decomposition, the most popular strategies are the Finite Volume Method (FVM) and the Discrete Ordinates Method (DOM) [START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF]. The angular decomposition is well suited to parallelism and obtains quite good efficiencies. However, it is limited to the number of discrete directions when the medium is gray or to the product of the number of discrete directions by the number of wavelengths when the medium is non-gray, which limits the degree of parallelism. Moreover, in such strategies the data for the whole spatial domain must be stored in every processor which may induce a memory limitation. An extensive review on the parallelization of the DOM and the FVM for solving the RTE is provided, for instance, by Coelho [START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF] (section 9

Parallel implementation). More recent works on the topic are given in [START_REF] Jolivet | Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning[END_REF][START_REF] Hu | The spatial and angular domain decomposition method for radiation heat transfer in 2D rectangular enclosures with discontinuous boundary conditions[END_REF][START_REF] Badri | High performance computation of radiative transfer equation using the finite element method[END_REF][START_REF] Lygidakis | Using a parallel spatial/angular agglomeration multigrid scheme to accelerate the FVM radiative heat transfer computation-Part I: Methodology[END_REF][START_REF] Lygidakis | Using a parallel spatial/angular agglomeration multigrid scheme to accelerate the FVM radiative heat transfer computation-Part II: Numerical results[END_REF][START_REF] Chen | Application of collocation spectral domain decomposition method to solve radiative heat transfer in 2D partitioned domains[END_REF]. We also found two older works on the subject [START_REF] Gao | Parallel multigrid solver of radiative transfer equation for photon transport via graphics processing unit[END_REF][START_REF] Yildiz | A parallel solution to the radiative transport in three-dimensional participating media[END_REF] which are not cited in [START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF]. Most parallel algorithms work either with MPI or openMP and sometimes with GPU [START_REF] Gao | Parallel multigrid solver of radiative transfer equation for photon transport via graphics processing unit[END_REF]. They use either angular decomposition or spatial decomposition. The spatial decomposition is commonly employed in CFD but its parallel efficiency is usually lower due to stronger dependencies between spatial elements. Recent developments have however investigated improved sweeping strategies, and the parallelization of alternative algorithms or alternative formulations [START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF][START_REF] Jolivet | Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning[END_REF][START_REF] Hu | The spatial and angular domain decomposition method for radiation heat transfer in 2D rectangular enclosures with discontinuous boundary conditions[END_REF][START_REF] Badri | High performance computation of radiative transfer equation using the finite element method[END_REF][START_REF] Lygidakis | Using a parallel spatial/angular agglomeration multigrid scheme to accelerate the FVM radiative heat transfer computation-Part I: Methodology[END_REF][START_REF] Lygidakis | Using a parallel spatial/angular agglomeration multigrid scheme to accelerate the FVM radiative heat transfer computation-Part II: Numerical results[END_REF][START_REF] Chen | Application of collocation spectral domain decomposition method to solve radiative heat transfer in 2D partitioned domains[END_REF]. By using a spatial decomposition of either Cartesian grids [START_REF] Yildiz | A parallel solution to the radiative transport in three-dimensional participating media[END_REF] or CFD unstructured grids [START_REF] Jolivet | Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning[END_REF][START_REF] Lygidakis | Using a parallel spatial/angular agglomeration multigrid scheme to accelerate the FVM radiative heat transfer computation-Part I: Methodology[END_REF][START_REF] Lygidakis | Using a parallel spatial/angular agglomeration multigrid scheme to accelerate the FVM radiative heat transfer computation-Part II: Numerical results[END_REF][START_REF] Nguyen | A massively parallel solution strategy for efficient thermal radiation simulation[END_REF], some numerical methods allow to reformulate the problem under a large linear system resulting from the discretized RTE [START_REF] Jolivet | Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning[END_REF][START_REF] Badri | High performance computation of radiative transfer equation using the finite element method[END_REF][START_REF] Chen | Application of collocation spectral domain decomposition method to solve radiative heat transfer in 2D partitioned domains[END_REF]. Then, some efficient parallel libraries such as Hypre or PETSc for example are used to solve the large linear system. As mentioned above, the spatial decomposition is very convenient when coupling the RTE with fluid dynamics solvers. In addition, multiresolution of either angular and/or spatial discretizations are sometimes used in multigrid approaches [START_REF] Jolivet | Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning[END_REF][START_REF] Lygidakis | Using a parallel spatial/angular agglomeration multigrid scheme to accelerate the FVM radiative heat transfer computation-Part I: Methodology[END_REF][START_REF] Lygidakis | Using a parallel spatial/angular agglomeration multigrid scheme to accelerate the FVM radiative heat transfer computation-Part II: Numerical results[END_REF][START_REF] Gao | Parallel multigrid solver of radiative transfer equation for photon transport via graphics processing unit[END_REF] in order to accelerate the convergence of the iterative process. However, our numerical method for solving the RTE is quite different from the previous ones as it consists in a modified FVM. Indeed, the radiative intensity at each grid node is explicitly computed for each discrete direction in function of the radiative intensities at neighboring upstream nodes for that direction. This formulation provides very accurate results in complex gaseous models. In this work, we carried out a parallelization of our algorithm for solving the RTE with a WSGG model, using MPI and openMP for coupling angular and spectral decompositions to obtain a parallel algorithm usable on current clusters. It is worth noticing that to the best of our knowledge, very few authors explicitly take into account the hierarchy of parallel contexts (distributed memory between nodes and shared memory inside nodes) present in current parallel systems, as in [START_REF] Poitou | Parallel computation for radiative heat transfer using the DOM in combustion applications: direction, frequency, subdomain decomposition, and hybrid methods[END_REF]. In this last work, the authors evaluate in particular a hybrid decomposition over angular, spectral and spatial domains. Although interesting speedups are obtained for larger numbers of processors, the efficiency is drastically degraded by the overhead due to the spatial decomposition that implies more iterations. However, the inclusion of a spatial decomposition in a parallel RTE solver is mandatory to achieve higher scalability, and is especially useful when the RTE is coupled with CFD. The inclusion of such a decomposition in our solver will be the subject of a future work.

Concerning the treatment of the spectral dependence of radiative properties of participating gases, in spite of their accuracy, line-by-line (LBL) calculations are too costly for engineering applications that involve coupled, multi-phenomena calculations. As such, there has been a continuous interest in the development of gas models that lead to computation efficiency and accurate spectral integration of the radiative transfer in problems with strong variations in the thermodynamic state. Among the non-gray gas models, the weighted-sum-of-gray-gases (WSGG) model [START_REF] Liu | The impact of radiative heat transfer in combustion processes and its modeling -with a focus on turbulent flames[END_REF][START_REF] Coelho | Radiative Transfer in Combustion Systems[END_REF][START_REF] Modest | Radiative Heat Transfer in Turbulent Combustion Systems: Theory and Applications[END_REF] is very likely the most widely used in engineering applications due to its relative simplicity and robustness. In past years, new developments of the model have relaxed some of its limitations, such as application to gas mixtures with varying molar concentration ratios and non-atmospheric pressures [START_REF] Bordbar | Pressure-dependent weighted-sum-of-gray-gases models for heterogeneous -mixtures at sub-and super-atmospheric pressure[END_REF] and non-gray boundaries [START_REF] Fonseca | Grouping the wall spectral bands for an effective computation of the radiative transfer in participating media bounded by non-gray walls[END_REF]. As shown in a recent study on coupled combustion calculations [START_REF] Consalvi | Assessment of engineering gas radiative property models in high pressure turbulent jet diffusion flames[END_REF], in comparison with more advanced global gas model, such as SLW/FSCK, the model can provide comparable accuracy, even if not being the most accurate, with less computation time.

Concerning the spatial integration of the RTE, the Monte Carlo method [START_REF] Liu | The impact of radiative heat transfer in combustion processes and its modeling -with a focus on turbulent flames[END_REF][START_REF] Coelho | Radiative Transfer in Combustion Systems[END_REF][START_REF] Modest | Radiative Heat Transfer in Turbulent Combustion Systems: Theory and Applications[END_REF] is often employed, but due to the huge CPU times, other methods are often preferred for combustion applications. The spherical harmonic method (SHM), the DOM and the FVM are among the most widely used at present [START_REF] Liu | The impact of radiative heat transfer in combustion processes and its modeling -with a focus on turbulent flames[END_REF][START_REF] Modest | Radiative Heat Transfer[END_REF][START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF]. They provide a good compromise between accuracy and computational requirements, and they are relatively easy to integrate in CFD codes. For combustion application, radiative heat transfer has to be simulated in 3D geometrically complex systems and the DOM / FVM are suitable since they are formulated for general control volumes. Many numerical methods have been developed for solving the RTE but, to the best of our knowledge, very few works have been devoted to the application of these methods to gas radiative transfer in complex geometries for real combustion application for which such methods are needed. A comprehensive review of the these numerical methods are available in Liu et al. [START_REF] Liu | The impact of radiative heat transfer in combustion processes and its modeling -with a focus on turbulent flames[END_REF] and Coelho [START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF].

The objective of the present study is to propose an accurate and efficient numerical method for solving the RTE in 3D complex geometry combined with a gas radiative property model for combustion applications. In the literature, the integration over the wavenumber spectrum is systematically done after solving the radiative intensities for each position r r r

and Ω Ω Ω direction [START_REF] Coelho | Radiative Transfer in Combustion Systems[END_REF][START_REF] Modest | Radiative Heat Transfer in Turbulent Combustion Systems: Theory and Applications[END_REF][START_REF] Modest | Radiative Heat Transfer[END_REF][START_REF] Coelho | Turbulence-radiation interaction: From theory to application in numerical simulations[END_REF][START_REF] Modest | The Treatment of Nongray Properties in Radiative Heat Transfer -From Past to Present[END_REF]. Thus, the RTE has to be solved for each wavenumber and leads, most of the time, to compute some geometrical quantities (related to the angular and spatial discretizations of the domain) for each wavenumber. In particular, the cost to compute these quantities can be high for an unstructured tetrahedral mesh. We call "Classical method" the one that leads to solving the RTE in this way. In this work, the calculations over the wavenumber spectrum and the angular space are reversed. This new way to solve the RTE is called "New method" and gives the same results as the "Classical method". This inversion of the computations hierarchy enables to save a part of the geometric computations related to each angular direction. So, this procedure presents the advantage of reducing the CPU time in both sequential and parallel contexts as it will be shown in the results section. To the authors' best knowledge, this comparison has not been carried before, so in this study this gain in CPU time will be evaluated.

The outline of the article is as follows. Details of the modified Discrete Ordinates Method (MDOM) applied to the RTE with a non-gray WSGG model and a description of its parallel implementation are given in section 2. The version of the MDOM employed here is an extension of the one presented in [START_REF] Trovalet | Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries[END_REF] where the finite volume method was used for spatial discretization. In [START_REF] Trovalet | Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries[END_REF], the MDOM was applied to a gray medium with an absorption coefficient constant in the medium.

In order to check the accuracy and efficiency of the numerical method, two test cases on cubic and cylindrical combustion chambers with non-isothermal and non-homogeneous gas mixtures are presented in section 3. The results are discussed and compared with those of other studies on similar topics. Also, performance and memory consumption of the parallel version are addressed in this section. Concluding remarks are finally offered in the final section.

2 Numerical treatment

Angular and spatial discretizations of the domain

The angular and spatial discretizations of the domain are carried out before solving the spectral RTE whatever the method used ("Classical" or "New"). The quadratures S n [START_REF] Modest | Radiative Heat Transfer[END_REF] and LDFE-SA-2 (denoted L 2 ) [START_REF] Jarrell | Discrete-ordinates quadrature sets based on linear discontinuous finite elements[END_REF] are used here for the discretization of the angular domain. L 2 provides higher accuracy than S 8 or S 12 . The total number of discrete angular directions is N dir . In this work, a grid structure usually used in CFD is adopted. As in [START_REF] Trovalet | Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries[END_REF], the computational spatial domain was divided into four-node tetrahedral elements using unstructured grids and a cell-vertex formulation was adopted for the control volumes. It consists in building control volumes around each node of the grid. Those are polyhedrons with quadrilateral faces. The V P control volume surrounding node P P P was built with all tetrahedrons having node P P P as vertex (we refer to [START_REF] Trovalet | Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries[END_REF] for the details of the construction of V P ). Γ P is the boundary of V P and is subdivided into N f surface panels. For a panel f , i f i f i f is the integration point located at the centre of gravity of the panel, A f is its surface area, and n f n f n f is the outward unit vector normal to the panel f .

DOM applied to the spectral RTE

By employing the non-gray WSGG model, the integration of the spectral RTE (for a non-scattering medium), over V P and within a ∆Ω k solid angle centred around Ω k Ω k Ω k direction yields [START_REF] Modest | Radiative Heat Transfer in Turbulent Combustion Systems: Theory and Applications[END_REF][START_REF] Modest | Radiative Heat Transfer[END_REF][START_REF] Trovalet | Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries[END_REF]:

N f f =1 A f I k ν,i f ∆ k f = κ ν,P a ν [T P ] I b [T P ] -I k ν,P ∆Ω k V P , (1)
where I k ν,P and I k ν,i f are approximations of the spectral radiative intensity in the Ω k Ω k Ω k direction at node P P P and the i f integration point, respectively.

∆ k f = Ω k Ω k Ω k ⋅ n f n f n f ω k
where ω k is the weight factor associated to Ω k Ω k Ω k according to the S n quadrature [START_REF] Modest | Radiative Heat Transfer[END_REF]. Here, κ ν,P and T P are respectively the spectral gray gas absorption coefficient and the local medium temperature at node P P P . The coefficient a ν [T P ] is a weight factor which depends on the local medium temperature at node

P P P . The function I b [T ]
is the blackbody intensity at the medium temperature.

Spatial discretization of the transport term of the spectral radiative intensity

To solve the set of equations ( 1), closure relations are needed between the I k ν,i f integration-point values and the nodal values of the spectral radiative intensity. The directional nature of radiation needed to be taken into account in order to establish the closure relations. Thus, for a specific direction of radiation, only the nodal values located upstream from the integration-point should be considered. It is known that the Step scheme causes strong false scattering while the Diamond scheme may yield oscillatory solutions [START_REF] Modest | Radiative Heat Transfer[END_REF][START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF]. The Exponential scheme relies on the integral form of the RTE and is expected to be more accurate than the two other schemes [START_REF] Modest | Radiative Heat Transfer[END_REF][START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF][START_REF] Trovalet | Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries[END_REF]. It has been then implemented in the present study.

The exponential scheme

With the Exponential scheme, the discrete equation ( 1) can be rewritten as [START_REF] Trovalet | Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries[END_REF]:

N f f =1 I k ν,u f D k ν,f = κ ν,P a ν [T P ] I b [T P ] -I k ν,P ∆Ω k V P - N f f =1 C k ν,f , (2)
with

D k ν,f = A f ∆ k f exp ⎛ ⎝ - i f u f κ ν (s) ds ⎞ ⎠ , (3)
C k ν,f = A f ∆ k f i f u f κ ν (s) a ν [T (s)] I b [T (s)] exp ⎛ ⎝ - i f s κ ν (u) du ⎞ ⎠ ds. (4)
The computation of ( 4) and (3) are presented in Appendix 1. It shoud be noticed that in the Step scheme, C k ν,f = 0 and andD k f for the Step scheme) are geometrical quantities that do not depend on the wavenumber. With the "New method", they are computed before carrying out the summation over the wavenumber spectrum.

D k ν,f = D k f = A f ∆ k f is independant of the wavenumber. Also A f ∆ k f and ∆Ω k V P (

Projections and linear interpolations

The projections and linear interpolations are taken from [START_REF] Trovalet | Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries[END_REF]. They were used to link I k ν,u f defined at the u f u f u f point with spectral radiative intensities at n b neighbouring upstream nodes from node P :

I k ν,u f ≅ 3 n b =1 α k f,n b I k ν,n b such that 3 n b =1 α k f,n b = 1 and 0 ≤ α k f,n b ≤ 1 (f = 1, 2, ..., 6 ; n b = 1, 2, 3). (5)
The way to compute the values for the coefficients α k f,n b is described in [START_REF] Trovalet | Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries[END_REF] and is not repeated here. This interpolation technique demonstrated to improve the closure relations and the accuracy of the results. It can be noticed that there are no projections with the Step scheme. The projections (in the case of the Exponetial scheme) and linear interpolations are geometric calculus that do not depend on the wavenumber. So, with the "New method", this procedure is carried out before the summation over the wavenumber spectrum.

Solution procedure 2.4.1 Marching order map

Once the angular discretization was fixed (according to the S n quadrature), a "marching order map" giving the optimal order in which the nodal nodes of the mesh should be visited was built for each discrete angular direction [START_REF] Trovalet | Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries[END_REF]. It consisted of numerically representing the direction of propagation of radiation (table of connectivity of radiation) through the computational spatial domain. It should be noticed that the direction of propagation is independent of the wavenumber.

Therefore, the marching map is built once for all whatever the method used ("Classical" or "New").

2.4.2 Expression of the spectral radiative intensity at P P P nodal node and in the Ω Ω Ω k direction Combining (2) and relation [START_REF] Coelho | Turbulence-radiation interaction: From theory to application in numerical simulations[END_REF] of Appendix 2, an algebraic equation is obtained relating I k ν,P to the spectral radiative intensities at neighbouring upstream nodes from node P P P :

I k ν,P = b k ν,P + g k ν,P d k ν,P , (6)
with

b k ν,P = - {n b ∈(Eu) k P } I k ν,n b {f ∈ψ k n b ,P } α k f,n b D k ν,f ; g k ν,P = κ ν,P a ν [T P ] I b [T P ] ∆Ω k V P - N f f =1 C k ν,f , (7)
d k ν,P = κ ν,P ∆Ω k V P + {f ∈Υ k P } α k f,P D k ν,f . (8)
In ( 6), I k ν,P is independent of downstream spectral radiative intensities and the method is fully explicit. Then, the total radiative intensities can be computed such as:

I k P = N i=1 I k i,P , I k κ,P = N i=1 κ i,P I k i,P . (9)
where the subscript i represents the i-th gray gas and N is the total number of gray gases for the considered WSGG model. The partial radiative heat flux incident to the medium boundary, the partial quantity G κ and the partial radiative heat flux in the Ω k Ω k Ω k direction are respectively defined by:

Q k inc,P = I k P n P n P n P ⋅ Ω k Ω k Ω k ω k if n P n P n P ⋅ Ω k Ω k Ω k > 0, G k κ,P = I k κ,P ω k , Q k r Q k r Q k r ,P = I k P Ω k Ω k Ω k ω k . (10)
Then, their total quantities at P P P nodal node can be easily obtained by summing [START_REF] Lygidakis | Using a parallel spatial/angular agglomeration multigrid scheme to accelerate the FVM radiative heat transfer computation-Part I: Methodology[END_REF] over the discrete angular directions.

The "New method" leads to compute these radiative quantities in this way.

Parallel computing and memory management

Work distribution In the initial version (Classical method), the wavenumbers were distributed over groups of the available machines. In each group, the angular directions were distributed over the machines and over the internal cores of each machine. As in the second test case (Cylindrical combustion chamber) presented in the results section, the work distribution of 25 gray gases and 48 directions (S 6 quadrature) over a system for example of 30 machines with 8 cores would lead to 5 groups of 6 machines. In each group, 5 gray gases are computed. And for each gray gas, the angular directions are distributed over the cores in the machines of the group. So, each machine would compute 5 gray gases and 8 directions. This parallel scheme works quite well but we have observed that a significant part of the computation time is spent in retrieving the geometrical information and performing the projections of integration points. Although some of those computations are distinct for each gray gas (wavenumber), a non negligible part is common to all wavenumbers. So, we have considered the possibility to take advantage of that by performing the common parts only once for the entire set of wavenumbers. To do that, a complete reorganization of the work distribution was necessary. Indeed, the geometrical computations involved are performed inside the solving of the RTE for each angular direction. So, the work distribution hierarchy must be inverted and the wavenumbers must be managed inside each angular direction as depicted in Fig. 1.

Finally, we obtain in Algorithm 1 a new work distribution in which the angular directions are distributed over the entire system (machines and cores inside machines) and the wavenumbers are managed inside each angular direction.

[Figure 1 In this parallel scheme, according to the size of the computing system, it would be possible to distribute the wavenumbers for one angular direction over a set of local cores of the machine, thus leading to nested parallelism (threads inside threads). Indeed, it depends on the number of machines, on the number of physical cores available in each machine as well as on the number of directions assigned to each machine. So, the local work distribution inside a single machine must be adapted according to the number of assigned angular directions and the total number of wavenumbers. At an extent, when the number of available machines is equal to the number of angular directions, they can be totally distributed over the machines (one direction per machine). Then, inside each machine only one thread is used to manage the unique direction and the remaining threads are used to compute the wavenumbers in parallel.

Let N loc dir be the number of angular directions assigned to a machine and N c the number of cores inside the machine.

Then, the number N td of threads used for the angular directions and the number N tw of threads used for the wavenumbers should satisfy the following conditions: N td × N tw = N c and N td should be a divider of N loc dir (strong constraints); N tw should be a divider of N (weak constraint). We differentiate strong and weak constraints according to their influence over the overall performance. As can be seen, it is important to use all the available cores of the machine and to get a good balancing of the angular directions over the threads as they are the largest tasks to process. As the computations related to a single wavenumber are much smaller, it is less problematic to have a slight imbalance of the wavenumbers over the threads. For the same reason, it is not interesting to distribute the wavenumbers over several threads when N is quite small.

In such cases, it is better to use all cores for the angular directions. Concerning the inter-machines data exchanges that take place at line 13 in Algorithm 1, they are limited to the radiative quantities

Q inc , G κ and Q r Q r Q r .
As these values are the results of integrations over the angular directions and the wavenumbers, the amount of data to transfer is only proportional to the mesh size and even less for Q inc (mesh surface). In fact, each machine computes partial integrations of these values according to the assigned angular directions. Then, in order to obtain the global results over all the machines, global reduction operations (additions) are necessary. These communication operations are performed only once per iteration of the RTE solving. So, the total amount of communications during the process stays quite limited.

The current design of our algorithm allows the use of at most N dir × N cores. This number stays quite limited when the number of gray gases is small although many parallel systems provide thousands of cores. So, it may be interesting to study the possibility to extend the parallelism of our algorithm. One solution would be to include a spatial domain decomposition method [START_REF] Poitou | Parallel computation for radiative heat transfer using the DOM in combustion applications: direction, frequency, subdomain decomposition, and hybrid methods[END_REF][START_REF] Plimpton | Parallel Sn Sweeps on Unstructured Grids: Algorithms for Prioritization, Grid Partitioning, and Cycle Detection[END_REF]. In our case, the spatial decomposition should be done for each angular direction. Nevertheless, the performance gain to expect is not totally obvious due to the sequential nature of the computations for each direction. Indeed, it would require additional computations (balanced decomposition for each angular direction) and some approximations in the computations to relax the dependency constraints between the vertices in the mesh in order to avoid frequent synchronizations between cores. However, such an extension would allow a better exploitation of large parallel systems and would deserve a complete study to determine its interest. As mentioned in the introduction, this will be the subject of a future work.

Memory management It must take into account the new work distribution described previously. In order to comply with this new distribution, several modifications have been done. First of all, data arrays have been re-organized in order to put the wavenumbers in the last dimension of the arrays. Indeed, as all wavenumbers are processed for each angular direction, it is better to store them in a compact way, in order to reduce cache misses at processor level. Moreover, in order to reduce the memory consumption on each machine, the sizes of the data arrays depending on the angular directions are limited to the number of active threads. This means that when the number N loc dir of angular directions assigned to a machine is larger than the number of threads N td , data will be stored only for N td directions at the same time. Then, the arrays are re-used for other angular directions through a sequential iteration other the angular directions. These modifications partly compensate the need to store data for all the wavenumbers. Indeed, the storage for each wavenumber must be distinct as the RTE must be solved independently for each of them. It is only once all wavenumbers are processed that integrations are done. With the WSGG model, this scheme is pertinent as the number of wavenumbers stays quite limited (up to 25).

Results and discussion

As stated before, the WSGG model was used in this work, considering several recent formulations namely gray and nongray models by Johansson et al. [START_REF] Johansson | Account for variations in the H2O to CO2 molar ratio when modeling gaseous radiative heat transfer with the weighted-sum-of-grey-gases model[END_REF], non-gray model by Dorigon et al. [START_REF] Dorigon | WSGG correlations based on HITEMP2010 for computation of thermal radiation in non-isothermal, non-homogeneous H2O/CO2 mixtures[END_REF] and non-gray model by Bordbar et al. [START_REF] Bordbar | A line by line based weighted sum of gray gases model for inhomogeneous CO2-H2O mixture in oxy-fired combustion[END_REF]. To determine the WSGG parameter sets, the HITRAN92 spectroscopic databases have been used in mixtures are presented hereafter for checking the implementations of these formulations. First, a gas mixture with uniform species concentrations (including air and oxyfuel combustions) in rectangular parallelepiped cavity is considered, whereas in the second test case a gas mixture with heterogeneous species concentration field in cylindrical combustion chamber is assumed. The comparison between the WSGG formulations of refs [START_REF] Johansson | Account for variations in the H2O to CO2 molar ratio when modeling gaseous radiative heat transfer with the weighted-sum-of-grey-gases model[END_REF], [START_REF] Dorigon | WSGG correlations based on HITEMP2010 for computation of thermal radiation in non-isothermal, non-homogeneous H2O/CO2 mixtures[END_REF], [START_REF] Bordbar | A line by line based weighted sum of gray gases model for inhomogeneous CO2-H2O mixture in oxy-fired combustion[END_REF] and the evaluation of these WSGG formulations against LBL and/or SNB calculations in 3D configurations has not been reported in the literature.

Rectangular parallelepiped cavity

The problem of air combustion where the gas is a mixture of 10 % CO2 and 20 % H2O was studied by Liu [START_REF] Liu | Numerical solutions of three-dimensional non-grey gas radiative transfer using the statistical narrow-band model[END_REF] and recently by Fraga et al. [START_REF] Fraga | Benchmark Solutions of Three-Dimensional Radiative Transfer in Nongray Media Using Line-by-Line Integration[END_REF]. The problem of oxyfuel combustion where the gas is a mixture of 85 % CO2 and 10 % H2O was studied by Porter et al. [START_REF] Porter | Evaluation of solution methods for radiative heat transfer in gaseous oxy-fuel combustion environments[END_REF]. For both air and oxyfuel combustions, the dimension of the cavity is 2 m × 2 m × 4 m. All the surrounding walls are black and cold at 300 K. The pressure of the gas mixture in the cavity is kept at 1 atm. The gas temperature is non-uniform but symmetrical about the centerline of the cavity and is specified in terms of T (r) = (T c -T e ) f (r) + T e where, r is the radial distance from the cavity centerline, T c is the gas temperature along the cavity centerline and T e is the exit temperature at z = 4 m. The variation of gas temperature is defined by f (r) = 1 -3 r 2 + 2 r 3 . The gas temperature outside the circular region is assumed to be uniform and at the value of the exit temperature. The centerline temperature T c is assumed to increase linearly from 400 K at the inlet z = 0 to 1800 K at z = 0.375 m, then decreases linearly to 800 K at the exit. For this simple 3D geometry, the spatial grid used by Liu The effect of angular discretisation on the results was performed using the S 6 and S 8 quadratures when the MDOM is used with the Exponential scheme and mesh 1. For the radiative heat source, it was found that the results obtained with S 6 and S 8 are very close and then only the results for S 8 are presented. The wall radiative heat flux along line 2 is slightly overpredicted with S 6 (in comparison with S 8 ). In the oxyfuel combustion case, the result for the wall radiative heat flux along line 3 is better with S 6 . Also, the innacurate representation of gas absorptivities by the WSGG model leads to more important errors in the radiative heat source in the oxyfuel combustion case.

Cylindrical combustion chamber

This second test case deals with a cylindrical combustion chamber studied numerically by Centeno et al. [START_REF] Centeno | Evaluation of gas radiation heat transfer in a 2D axisymmetric geometry using the line-by-line integration and WSGG models[END_REF]. A natural gas combustion is used in a 600 kW furnace [START_REF] Garréton | First Aerodynamics of Steady State Combustion Chambers and Furnaces Workshop[END_REF][START_REF] Centeno | The influence of gas radiation on the thermal behavior of a 2D axisymmetric turbulent non-premixed methane-air flame[END_REF] to generate a realistic non-isothermal and non-homogeneous flame.

The measured gas temperature profiles and chemical species concentrations were presented in [START_REF] Garréton | First Aerodynamics of Steady State Combustion Chambers and Furnaces Workshop[END_REF], which had been used by other researchers [START_REF] Centeno | Evaluation of gas radiation heat transfer in a 2D axisymmetric geometry using the line-by-line integration and WSGG models[END_REF][START_REF] Centeno | The influence of gas radiation on the thermal behavior of a 2D axisymmetric turbulent non-premixed methane-air flame[END_REF][START_REF] Magel | Simulation of detailed chemistry in a turbulent combustor flow[END_REF]. The dimensions of the chamber are representative to those found in industrial systems.

The cylinder has a radius r of 0.25 m and a length of 1.7 m (Fig. 9(a)).

[Figure 9

The coefficients for the correlation [START_REF] Lygidakis | Using a parallel spatial/angular agglomeration multigrid scheme to accelerate the FVM radiative heat transfer computation-Part II: Numerical results[END_REF] are given in Table 1. The correlation ( 11) is valid for the whole chamber domain, and provides the value of the mean temperature (in K) are every position (r and z should be in units of meter).

The H2O and CO2 molar distributions were fitted by trigonometric correlations, each one valid for the whole chamber domain, as given by the relations below [START_REF] Centeno | Evaluation of gas radiation heat transfer in a 2D axisymmetric geometry using the line-by-line integration and WSGG models[END_REF]:

X H2O (z, r) = a w cosh ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ π (r -b w ) c w ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ sin ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ π (z -d w ) f w ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ + g w , (12)
X CO2 (z, r) = a c cosh ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ π (r -b c ) c c ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ sin ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ π (z -d c ) f c ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ + g c , (13)
The coefficients for these correlations are also given in Table 1. As for the temperature distribution, r and z should be in units of meter. The molar distributions ( 12) and ( 13) are used in the considered WSGG model for the calculation of the gas absorption coefficient and the weight factor.

[Table 1 about here.]

The present results are compared with those obtained by Centeno et al. [START_REF] Centeno | Evaluation of gas radiation heat transfer in a 2D axisymmetric geometry using the line-by-line integration and WSGG models[END_REF]. Profiles of the radiative heat source along the chamber centerline and the radial direction at the axial positions z = 1 m and z = 1.5 m are displayed in Fig. 10. The reference solution was obtained with two WSGG models to compute radiative heat transfer in the H2O and CO2 mixtures as well as the LBL calculations (for comparison/validation). The first WSGG model named CR-WSGG considers correlations for a constant ratio between the H2O and CO2 molar concentrations (X H2O X CO2 = 2), while the second named NCR-WSGG model makes use of separate correlations for H2O and CO2, which are combined for any concentration ratio. The WSGG correlations as well as the LBL calculations were based on the HITEMP2010 database. The number of gray gases, for the CR-WSGG and NCR-WSGG models, is equal to 5 and 25 respectively. As the LBL calculations are very time consuming for 3D geometries and for the study of coupled heat transfers, only the CR-WSGG and NCR-WSGG models were implemented in this work. The WSGG coefficients are given in Centeno et al. [START_REF] Centeno | Evaluation of gas radiation heat transfer in a 2D axisymmetric geometry using the line-by-line integration and WSGG models[END_REF] and are not repeated here.

In Centeno et al. [START_REF] Centeno | Evaluation of gas radiation heat transfer in a 2D axisymmetric geometry using the line-by-line integration and WSGG models[END_REF], the spatial discretisation of the RTE was accomplished with the DOM using the Diamond scheme to compute intensities in terms of upstream values. It is well known that the DOM leads to some errors in the numerical Our results are however in good agreement with those of Centeno et al. For the same WSGG model and with the finer spatial grid used in this case, very few discrepencies can be seen between the results obtained with the Exponential and

Step schemes for the radiative heat source along the chamber centerline (Fig. 10(a)) and the radial direction at position z = 1.5 m (Fig. 10(c)). We carried out other simulations (not presented here) and it was found that notable discreapencies were obtained when the spatial grid is coarser. As it is expected, the result are more accurate with the NCR-WSGG model in comparison with the CR-WSGG model.

[Figure 10 about here.]

The numerical simulations were performed with a 16 cores computer (2 CPUs Intel @2.0GHz,8 cores/CPU, 252GBRAM) using Hyper-Threading and Intel C compiler. Table 2 presents the CPU times between the different numerical methods on the same mesh composed of 570,456 nodes (3,361,362 tetrahedrons) and the S 6 angular discretization. In this Table , only the angular directions are computed in parallel using 16 threads. As expected, the Step scheme (in comparison with the Exponential scheme) and the CR-WSGG model (in comparison with the NCR-WSGG model) are less time consuming. In both cases the transparent windows are being accounted as one of the gray gases (with absorption coefficient equal to zero) since the temperatures on the boundaries are not uniform. It can be seen that the so called "New method" is faster than the "Classical method". Moreover, as the parallelism is only applied to process 48 directions with 16 threads (3 directions per thread), this means that the ratios in Table 2 are directly linked to the maximal CPU times to process three directions in sequential. Also, as the computation amounts are very similar between distinct directions, the CPU times between the threads are very closed to each other. This implies that the ratios between those parallel executions are representative of the ratios that would be obtained in purely sequential executions. So, those results also show that the "New method" not only improves the CPU time of our RTE solver in parallel executions but also in purely sequential ones. The ratio is always higher when the Step scheme is used. The CPU time is significantly improved for the NCR-WSGG model where the ratio is equal to 16.18 in this case. The ratio increases with the number of gray gases used in the WSGG model: 16.18 againts 3.64 (Step scheme) and 2.85 againts 1.79 (Exponential scheme). It is also interesting to see that with the New method and the Step scheme, the computation time increases very slightly with the number of gray gases: 39 s. for 5 gray gases (CR-WSGG model) against 43 s. for 25 gray gases (NCR-WSGG model).

[Table 2 about here.]

Performance evaluation of parallelism and memory consumption

The overall performance of parallelism as well as the memory consumption have been evaluated for the second test case (Cylindrical combustion chamber) and with the exponential scheme. For this part, all the simulations have been done on the Grid'5000 plateform [START_REF] Cappello | Grid'5000: a large scale, reconfigurable, controlable and monitorable Grid platform[END_REF][START_REF]Grid'5000 platform[END_REF]

Performance evaluation of parallelism

As there are two levels of parallelism, it was pertinent to check how efficient is each of them. So, we have performed three series of simulations. The first one uses a single machine and an increasing number of threads. The second one uses only one thread per machine and an increasing number of machines. The last one is a combination of machines and threads in order to maintain the total number of used cores equal to the number of discrete angular directions (here equal to 48). As can be seen in Fig. 11(a), the multi-threading stays quite close to optimal efficiency up to 24 threads. Although the machine contains 32 physical cores, the number of cores used in this simulation is limited to 24 as it is the largest number under 32 that divides the number of discrete angular directions (48). So, using 32 threads in this case would not bring any benefit over the performance as some of the cores would still have to compute two directions one after the other.

The slight decrease observed when the number of cores increases is quite a classical phenomenon. It comes from the fact that with more threads, the overall computation time decreases while the parts of the program that are not subject to parallelism still take the same time. So, these sequential parts of the program become more significant and prevent to get an exact division of the computation time by the number of threads.

In Fig. 11(b), both speedup and efficiency stay close to the optimal value when the number of machines increases. This is a very good result as it shows that the communications induced by the data distribution do not imply any significant performance loss. This comes from the fact that the only data exchanged between the machines are the partially integrated values of the heat fluxes. So, their respective amounts are limited to the mesh size, which is quite small compared to the overall amount of data managed inside each machine.

Finally, in Fig. 11(c) are given the computation times for a combination of multi-machine and multi-thread parallelism. In this context, the 48 angular directions are distributed over the machines and inside each machine, the appropriate number of threads is used to manage the local directions in parallel. For example, with 3 machines, each machine is in charge of 16 directions and must then use 16 threads. The interest of this simulation is to determine whether there is a type of parallelism to choose in priority. The results indicate that the performance is slightly better when using more machines and less threads. This result is quite counter-intuitive as generally, performance is mainly impacted by the communications.

However, in our case, we have already seen that the communications have little impact over the performance. The main degradation comes from the increase of the number of threads. This may be due to several factors related to the memory usage. As the threads share the same memory space, some particular memory accesses may slow down their executions despite the care taken in the way the arrays are built. Also, when the number of threads increases, the amount of memory used also increases and this may induce some slow down in the memory accesses.

Memory consumption

In Fig. 11(d) are depicted the memory consumptions in the three parallel situations described previously. In the context of a single machine (SM-MT), the memory consumption linearly increases with the number of threads. This is expected as in order to ensure the coherency of the computations in the shared memory between threads, they must have some private copies of a part of the data arrays. However, those copies represent a quite small part of the entire data set and this can be seen in the slope of the curve that is much smaller than 1. Indeed, the memory consumption is multiplied only by 2 when the number of threads is multiplied by 12.

The last two cases correspond to the situations with multiple machines. Using a single thread per machine (MM-ST) is the most economical strategy in terms of memory consumption. Indeed, in this case there is no array duplication due to several threads. Moreover, the memory requirements decrease when the number of machines increases due to the distribution of the angular directions over the machines. This phenomenon can also be observed in a stronger way in the last parallel situation in which multiple threads are used inside each machine (MM-MT). In fact, with 2 machines, there are 24 threads per machine and the memory amount is thus merely the same as in the SM-MT situation with 24 threads.

Then, when the number of machines increases, the number of threads per machine decreases producing a sharp memory saving as there are less array duplications per machine.

Finally, we obtain a multi-level parallel algorithm that can efficiently takes advantage of multi-scale parallel systems.

Also, the special care brought to memory consumption provides the ability to treat large size problems on parallel systems with a limited number of machines.

[Figure 11 about here.]

Conclusion

This paper reports a new numerical solution of radiative heat transfer in a 3D combustion system with a gas mixture using a parallel modified discrete ordinates method (MDOM) applied to 3D unstructured tetrahedral meshes. Several recent formulations of the weighted-sum-of-gray-gases (WSGG) model were implemented. To reduce the false scattering, the MDOM was combined with a closure relation based on an exponential scheme and was compared to the Step scheme. Our results were then compared with published results on two cases with non-isothermal gas mixtures. First, a gas mixture with uniform species concentrations (including air and oxyfuel combustions) was considered, whereas in the second test a gas mixture with heterogeneous species concentration field was assumed. The use of non-gray WSGG models showed a very good level of accuracy for the evaluation of wall heat transfer and radiative heat source. The Exponential scheme is more accurate than the Step scheme, in particular when the spatial grid is coarser. The summation over the wavenumber spectrum of the radiative model is usually carried out after computing the radiative quantities for each direction. By summing first over the wavenumber spectrum, a reduction in the computation time was obtained. In the second test case, the acceleration factor is 16 for the Step scheme which gives very close results with the Exponential scheme. Moreover, the use of a hierarchical parallel algorithm combining a work distribution over machines and multiple threads inside each machine allows us to obtain very good performance for the most demanding test case (cylindrical combustion chamber).

This demonstrates the efficiency of the implemented parallel scheme. With the New method and the Step scheme, the computation time increases very slightly with the number of gray gases: 39 s. for 5 gray gases (CR-WSGG model) against 43 s. for 25 gray gases (NCR-WSGG model). It should be noted that the procedure to speed up the calculations is not specific to our method. It can be applied to any numerical method for solving the spectral RTE.

In future research, we plan to study the turbulence and radiation interactions (TRI). Radiative heat transfer appears then as a source term in the overall energy conservation equation and is only one of the many complex phenomena that must be modeled in TRI. The present numerical method and the so called "LS-STAG code" [START_REF] Nikfarjam | The LS-STAG immersed boundary/cut-cell method for non-Newtonian flows in 3D extruded geometries[END_REF][START_REF] Portelenelle | Accurate discretization of diffusion in the LS-STAG cut-cell method using diamond cell techniques[END_REF] will be used as radiative and Indeed, if f ∈ Ψ k P then f ∈ ψ k n b ,P if n b ∈ E k f,P . Then, according to [START_REF] Trovalet | Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries[END_REF], the left member of (2) can be expressed as: 11), and chemical species distributions, Eqs. ( 12) and ( 13), as given in Centeno et al. [START_REF] Centeno | Evaluation of gas radiation heat transfer in a 2D axisymmetric geometry using the line-by-line integration and WSGG models[END_REF].

N f f =1 I k ν,u f D k ν,

A surface area, m 2 a 1 N

 21 ν weight factor in the non-gray WSGG model G G G isobarycentre i i i integration point on surface panel I total radiative intensity, W m -2 sr -number of gray gases N dir number of discrete angular directions N f number of surface panels of i f integration point in an upstream position

2. 4 . 3

 43 Expression of heat fluxes at P P P nodal node and in the Ω Ω Ω k direction

  about here.] Algorithm 1: Global parallel scheme Distribute the angular directions over the machines foreach machine do in parallel foreach thread inside current machine do in parallel Choose one Ω Ω Ω angular direction among the angular directions assigned to current machine foreach P P P nodal node in the mesh along Ω Ω Ω direction do foreach wavenumber ν do (possibly in parallel) Solve the RTE for P P P , Ω Ω Ω and ν done done done Compute the partial integrations of the radiative quantities related to the assigned directions Reduce the partial integrations by summation between machines (communication) done done

  and also by Porter et al. is structured and is 17 × 17 × 24. It is uniform in the x and y directions and non-uniform in the z direction with finer grids placed around the peak gas temperature (z = 0.375 m). It is referred as mesh 1. Fraga et al. used a finer structured spatial grid composed of 64 × 64 × 128 nodes and is uniform in the three directions x, y and z. It is referred as mesh 2. We used both the spatial meshes in our simulations. Benchmark data of Liu and Porter et al. were obtained using the ray tracing method with T 7 quadrature and the Statistical Narrow-Band (SNB) model. Those of Fraga et al. were obtained using a FVM with high angular resolution and the LBL model with the more up-to-date spectroscopic database HITEMP 2010. In the present work, the MDOM using the Step and Exponential schemes were used. The results for the S 6 (48 directions), S 8 (80 directions) and L 2 (128 directions) quadratures are presented. The radiative heat source along the centerline (1 m, 1 m, z) named line 1 and the wall radiative heat flux along line 2 (2 m, 1 m, z) are shown in Figs. 2, 3, 4 for the air combustion case and in Figs. 5, 6, 8 for the oxyfuel combustion case. The wall radiative heat flux along line 3 (x, 1 m, 4 m) is only shown in Fig. 7 for the oxyfuel combustion case since no benchmark data is available for the air combustion case. In Figs. 2 and 3, differences are seen between the Liu's data with those of of Fraga et al. They are more pronounced for the wall radiative heat flux since it is more sensitive to the calculations. This is because the solution computed by Fraga et al. was based on LBL integration and HITEMP 2010 database (the same database used to generate the coefficients in the non-gray model by Dorigon et al.) whereas Liu used Narrow-Band model based on a different database. The data by Fraga et al. are then more acccurate and they are considered as the reference.The comparison of our results with the reference solution for the air combustion case is shown in Figs.2, 3, 4, which includes the results for the gray and non-gray models by Johansson et al.[START_REF] Johansson | Account for variations in the H2O to CO2 molar ratio when modeling gaseous radiative heat transfer with the weighted-sum-of-grey-gases model[END_REF] and the non-gray model by Dorigan et al.[START_REF] Dorigon | WSGG correlations based on HITEMP2010 for computation of thermal radiation in non-isothermal, non-homogeneous H2O/CO2 mixtures[END_REF]. As expected, the non-gray models are in better agreement with the reference solution. Similar conclusions can be drawn for oxyfuel combustion (Figs.[START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media[END_REF][START_REF] Novo | Parallelization of the Discrete Transfer Method[END_REF][START_REF] Jolivet | Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning[END_REF][START_REF] Hu | The spatial and angular domain decomposition method for radiation heat transfer in 2D rectangular enclosures with discontinuous boundary conditions[END_REF] where, the results using the gray and non-gray models by Johansson et al. and the non-gray model by Bordbar et al.[START_REF] Bordbar | A line by line based weighted sum of gray gases model for inhomogeneous CO2-H2O mixture in oxy-fired combustion[END_REF] together with the reference solution are presented.For the two combustion cases (air and oxyfuel), the results obtained with the non-gray WSGG model are the closest to the reference solution when the MDOM is used with the Exponential scheme. Although the MDOM with the Step scheme has the advantage to be faster (about three CPU times faster than the MDOM based on the Exponential scheme), it tends to underpredict the radiative heat source (Figs. 2(a), 5(a)) and the wall radiative heat flux (Figs. 3(a), 6(a), 7(a)) leading to important errors.

Fig. 4 (

 4 Fig. 4 (air combustion case) shows that our results obtained with the non-gray models by Johansson et al. and Dorigon et

Fig. 8 (

 8 Fig. 8 (oxyfuel combustion case) shows our results obtained with the non-gray models by Johansson et al. and Bordbar

  about here.] The fuel composed of 90 % of methane and 10 % of nitrogen (by mass fraction) is injected through a central cylindrical tube having 0.03 m radius. The mass flow rate of fuel is 0.01453 kg/s at a temperature of 313.15 K. The air composed of 23 % of oxygen, 76 % of nitrogen and 1 % of water vapor is supplied by an annular duct having 0.05 m outer radius. The mass flow rate of air is 0.1988 kg/s at a temperature of 323.15 K. The temperature of the wall is fixed at 393.15 K, and the operating pressure is 101325 Pa.The cylinder side wall and ends are assumed as black surfaces and are kept at a temperature of 393.15 K. The two ends of the geometry are inlet and exit openings for air, fuel, and flue gases (Fig.9(a)). The experimental temperature field was fitted by a fifth degree polynomial function correlation as given by[START_REF] Centeno | Evaluation of gas radiation heat transfer in a 2D axisymmetric geometry using the line-by-line integration and WSGG models[END_REF]:T (z, r) = c 1 + c 2 r + c 3 z + c 4r 2 + c 5 r z + c 6 z 2 + c 7 r 3 + c 8 r 2 z + c 9 r z 2 + c 10 z 3 + c 11 r 4 + c 12 r 3 z +c 13 r 2 z 2 + c 14 r z 3 + c 15 z 4 + c 16 r 5 + c 17 r 4 z + c 18 r 3 z 2 + c 19 r 2 z 3 + c 20 r z 4 + c 21 z 5 ,

simulations. The S 6 (

 6 48 directions) was employed for angular discretization of the RTE both in Centeno et al. and in this our simulations. A 2D axisymmetric grid was used consisting of 50 and 90 cells in the radial and axial directions. It should be noted that the differences between our results and those of Centeno are due to the model and the spatial discretisation of the RTE. We used a real 3D model whereas a 2D axisymmetric model was used in Centeno et al. In Centeno et al, the cylinder side wall and ends are kept at a temperature of 393.15 K. In our case, they are kept at the wall temperature where the value are obtained by the polynomial function given in Eq. (11) for the corresponding wall position leading to closer results with Centeno et al. This behavior was also observed in the PhD thesis of Miranda [35]. In Centeno et al., the spatial discretisation of the RTE was accomplished with the DOM using the Diamond scheme whereas we used a MDOM with the Exponential and Step schemes. Finally, an unstructured tetrahedral spatial mesh composed of 570,456 nodes (3,361,362 tetrahedrons) was used to generate mesh independant results (Fig. 9(b)).

Figs. 11 (

 11 Figs. 11(a) and 11(b) show the speedup and the efficiency of the parallelism. The speedup indicates how much faster

Figure 1 :

 1 Figure 1: Work distribution over machines and threads.

Figure 4 : 8 Figure 5 :

 485 Figure 4: Air combustion. Comparison to benchmark data of radiative heat source along line 1 and wall radiative heat flux along line 2. The WSGG models of Non-gray Johansson et al. and Dorigon et al. were implemented. The MDOM was used with the Exponential scheme and a finer grid both in space (64 × 64 × 128 nodes) and angle (L 2 -128 discrete directions).

Figure 6 :

 6 Figure 6: Oxyfuel combustion. Comparison to benchmark data of wall radiative heat flux along line 2. Several methods for the WSGG model were implemented. The MDOM was used with the Step scheme, the Exponential scheme, S 6 and S 8 . 31

  (a) axial direction at centerline (b) radial direction at z = 1 m (c) radial direction at z = 1.5 m

Figure 10 :

 10 Figure 10: Comparison to benchmark data of radiative heat source along the chamber centerline and along the radial direction at the axial positions z = 1 m and z = 1.5 m. The CR and NCR WSGG models were implemented. The MDOM was used with the Exponential and Step schemes.

  (a) Speedup and efficiency in function of the number of threads on a single machine. (b) Speedup and efficiency in function of the number of machines with a single thread. (c) Time ratio for different combinations of machines and threads providing 48 cores. (d) Memory consumption in multi-threading, multi-machine and parallel combination contexts.

Figure 11 :

 11 Figure 11: Performance and memory consumption of the parallel scheme in different contexts.

  c

  Johansson et al., whereas the HITEMP 2010 databases have been used in Dorigon et al. and Bordbar et al. Two test cases with non-isothermal gas

  with machines under Linux Debian 4.19.181-1 containing 2 CPUs Intel Xeon Gold 6130 with 16 cores/CPU, 192GB RAM and 10Gb Ethernet. The C compiler used is gcc (Debian 8.3.0-6). All results are means of a set of executions.

Table 1 :

 1 10 0 2.7637 10 0 2.0665 10 -1 1.743 10 -22.4736 10 -1 3.0427 10 -1 d c f c g c 1.5762 10 0 2.6320 10 0 9.6858 10 -2 Coefficients for the temperature, Eq. (

	1	c 2	c 3	c 4	c 5	c 6
	301.9 10 0	-1.032 10 4	1228 10 0	4.467 10 5	2215 10 0	-5095 10 0
	c 7	c 8	c 9	c 10	c 11	c 12
	-4.409 10 6 -7.247 10 4	2.442 10 4	6286 10 0	1.724 10 7	6.402 10 5
	c 13	c 14	c 15	c 16	c 17	c 18
	-5.288 10 4 -2.765 10 4	-2011 10 0	-2.404 10 7	-1.493 10 6	5.978 10 4
	c 19 2.257 10 4	c 20 7156 10 0	c 21 -4.257 10 0 3.6337 10 -2 2.3154 10 -1 2.8847 10 -1 a w b w c w
	d w	f w	g w	a c	b c	c c
	1.5974

(a) MDOM using the Step scheme and S 6 (b) MDOM using the Exponential scheme and S 8

(a) Radiative heat source along line 1 (b) Wall radiative heat flux along line 2

Figure 2: Air combustion. Comparison to benchmark data of radiative heat source along line 1. Several methods for the WSGG model were implemented. The MDOM was used with the Step scheme, the Exponential scheme, S 6 and S 8 .

Appendix 1. Calculations of the source term in the RTE

In this section, the wavenumber is omitted to simplify the writing. For s ∈ (u f , i f ) and u ∈ (s, i f ), the absorption coefficient is evaluated as κ(s) ≈ κ = 1 2 κ i f + κ u f and κ(u) ≈ κ. Then,

where Λ f = i fuf is the distance between points i f and u f .

Let be f (s) = a[T (s)] I b [T (s)] for s ∈ (u f , i f ). We have to compute the following integral:

In that case, the function f is assumed to vary linearly along the optical path (u f , i f ):

Then, it follows:

a[T i f ] is computed from temperatures of the four vertices of the tetrahedron and the α i (i = 1, 2, 3, 4) coordinates of the i f integration points. Then, it can be:

Appendix 2. Expression of the left member of ( 2)

Let E k f,P be the set of vertices of a T tetrahedron located upstream from the i f integration point (on surface panel f ), such that i f ∈ Γ P and i f ∈ T . Let Υ k P be the set of surface panels f such that P belongs to E k f,P and let Ψ k P be the set of surface panels f such that P doesn't belong to E k f,P . If (E u ) k P is defined as the set of neighbouring upstream nodes from node P , then we can write:

Let,

,n b by the closure relations seen previsouly. It should be noted that ψ k n b ,P can also be defined from E k f,P and (E u ) k P . List of Tables

List of Figures

2.85

Table 2: CPU time and ratio between the so called "Classical method" and "New method".