Lorella Palluotto
email: lorella.palluotto@centralesupelec.fr

Nicolas Dumont

Pedro Rodrigues

Olivier Gicquel

Ronan Vicquelin

Assessment of randomized Quasi-Monte Carlo method efficiency in radiative heat transfer simulations

Keywords: Monte Carlo, Quasi-Monte Carlo, low-discrepancy sequence, radiative heat transfer, combustion

Radiation can play a central role in turbulent reactive flows where heat transfer is enhanced in applications with high temperature and pressure. The Monte Carlo method is a successful technique to solve the radiative transfer equation accurately with relative ease while retaining detailed properties.

However, its drawback is associated to a slow convergence rate. One strategy to improve the efficiency of Monte Carlo method consists in replacing the pseudo-random sequences with an alternative sampling: the low-discrepancy sequences. The introduction of such sequences in Monte Carlo leads to Quasi-Monte Carlo methods. Their advantage lies in a higher convergence rate compared to MC methods which have however not been assessed in 3D participating media. Additionally, in order to get an error estimation which is necessary in practical applications, a randomization of Quasi-Monte Carlo is needed (Randomized-QMC). Such Randomized-QMC methods have not been considered for simulations of radiative heat transfer in participating media before. In the present study, Monte Carlo and Randomized Quasi-Monte Carlo methods are assessed in terms of efficiency and computational cost in radiative heat transfer simulations of three practical 3D configurations. Comparisons in terms of local standard deviation, convergence rate, and final computational cost show that Randomized Quasi-Monte Carlo outperforms Monte Carlo in all the investigated cases.

Introduction

In turbulent flames, thermal radiation, turbulence and kinetics are fully coupled to each other in a highly non-linear way and the incorporation of all these contributions in numerical simulations is therefore a great challenge. In combustion applications, radiation can play a significant role, since the heat transfer from the flame to the walls is driven, apart from convection, also by radiation of burnt gases inside the chamber. Therefore, coupled simulations involving combustion and radiative heat transfer are more and more used and targeted. In order to correctly account for both convective and radiative contributions to wall fluxes, the simultaneous solution of the radiative transfer equation (RTE) and the governing equations for reactive flows is required. However accurate modeling of radiative heat transfer requires instantaneous spatially detailed information about temperature, pressure and species concentration. On the one hand, Reynolds Average Navier-Stokes (RANS) simulations does not provide directly such information since only spatially averaged results are available, with a relatively low computational cost. On the other hand, Direct Numerical Simulations (DNS) fully resolve the flow field in time and space, but such an approach remains unaffordable in most applications. In between, there are the Large Eddy Simulations (LES) that solve in time and space the energetic turbulent eddies while only modelling the small eddies. Coupled simulations based on LES are therefore a promising approach to describe accurately radiative heat transfer in turbulent flows.

Several families of methods exist to solve numerically radiative transfer [START_REF] Chandrasekhar | Radiative transfer[END_REF][START_REF] Modest | Radiative heat transfer[END_REF] such as Discrete Ordinates Method (DOM), the Finite Volume Method (FVM) and the spherical harmonics method. Another kind of method, which is the one considered in this study, is Monte Carlo (MC) integration, which is strongly appreciated for its ability to provide accurate results that can serve as a reference to others methods. However, its main drawback is the associated computational cost which can become prohibitive.

In coupled combustion simulations, different approaches are then combined to model the reactive turbulent flow and thermal radiation all together. First studies on radiation-combustion coupling were carried out in RANS simulations thanks to the low computational cost of such an approach. RANS simulations have been coupled to DOM for simulations of industrial-scale furnaces in [START_REF] Adams | Modeling effects of soot and turbulence-radiation coupling on radiative transfer in turbulent gaseous combustion[END_REF] with the gray gas approximation and [START_REF] Habibi | Impact of radiation models in cfd simulations of steam cracking furnaces[END_REF] using the Weighted Sum of Gray Gases (WSGG) model for radiative properties of gases. Coupled simulations of turbulent flames of methane and air can be found in [START_REF] Coelho | Spectral radiative effects and turbulence/radiation interaction in a non-luminous turbulent jet diffusion flame[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], using the WSGG approach. RANS computations are coupled to MC method in simulation of turbulent flames in [START_REF] Snegirev | Statistical modeling of thermal radiation transfer in buoyant turbulent diffusion flames[END_REF] with the help of global models for participating gases. More accurate models for radiative properties of gases are used in coupled RANS-MC simulations of turbulent flames in [START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF][START_REF] Zhao | A transported probability density function/photon monte carlo method for high-temperature oxy-natural gas combustion with spectral gas and wall radiation[END_REF][START_REF] Ren | Monte carlo simulation for radiative transfer in a high-pressure industrial gas turbine combustion chamber[END_REF] where computational domains of a few thousands cells are retained. Large-eddy simulations for turbulent flows have been coupled to DOM method in simulation of combustion chambers in [START_REF] Jones | Combination of dom with les in a gas turbine combustor[END_REF] using a gray gas model and in [START_REF] Poitou | Analysis of the interaction between turbulent combustion and thermal radiation using unsteady coupled les/dom simulations[END_REF][START_REF] Berger | On the sensitivity of a helicopter combustor wall temperature to convective and radiative thermal loads[END_REF] using global models for gases. The effects of subgrid-scale turbulence-radiation interactions have recently been investigated with LES coupled to an FVM solver of thermal radiation based on WSGG.

In some recent works, thanks to the increase in computing power and advances in numerical algorithms and solver scalability, coupled LES or DNS with MC methods have been achieved to carry out state-of-the-art simulations of radiative heat transfer. The applications have grown from the study of a statistically 1D turbulent premixed flame with a grey-gas assumption [START_REF] Wu | Direct numerical simulation of turbulence/radiation interaction in premixed combustion systems[END_REF] to the accounting of accurate gas radiative properties in a turbulent jet flame [START_REF] Gupta | Turbulence-radiation interactions in large-eddy simulations of luminous and nonluminous nonpremixed flames[END_REF], in turbulent boundary layers [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF] and a confined turbulent flame representative of an actual combustor [START_REF] Koren | Multiphysics simulation combining large-eddy simulation, wall heat conduction and radiative energy transfer to predict wall temperature induced by a confined premixed swirling flame[END_REF]. The use of MC methods in 3D unsteady simulations such as DNS and LES has therefore become affordable and enables high-fidelity of the described radiative heat transfer. However, such coupled and multiphysics simulations remain very costly and additional efforts are necessary to make MC methods more efficient.

The purpose of this study is then to investigate improvements of MC methods while aiming at unsteady scale-resolving simulations of turbulent combustion applications. Thanks to their statistical estimation of the radiative power and associated error control, Monte Carlo methods provide reference results. Moreover, the accounting for spectral gas radiative properties and complex geometries is straightforward. The drawback of all these methods is nonetheless their slow convergence, making them very much computationally intensive. Several improvements of physical and mathematical natures have been developed in the literature. For example, embedding the reciprocity principle into the Monte-Carlo simulation as in [START_REF] Walters | Rigorous development for radiation heat transfer in nonhomogeneous absorbing, emitting and scattering media[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] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF] reduces the conventional Monte Carlo convergence time compared to the Forward (FM) and Backward (BM) MC in the treatment of near-isothermal systems. From a mathematical point of view, variance reduction techniques [START_REF] Delatorre | Monte carlo advances and concentrated solar applications[END_REF][START_REF] Lemieux | Monte carlo and quasi-monte carlo sampling[END_REF] such as importance sampling have been considered in computations of radiative heat transfer [START_REF] De Lataillade | A net-exchange monte carlo approach to radiation in optically thick systems[END_REF][START_REF] Juvela | Efficient monte carlo methods for continuum radiative transfer[END_REF][START_REF] Feldick | Variance reduction techniques for monte carlo solution of radiative transfer in hypersonic flows: Hybrid p-1-monte carlo[END_REF][START_REF] Zhang | Optimized emission-based reciprocity monte carlo method to speed up computation in complex systems[END_REF] to strongly accelerate MC convergence by favouring more meaningful photon rays in the computed quantity of interest. An additional improvement is the use of an alternative sampling mechanism for numerical integration usually referred to as Quasi-Monte Carlo (QMC) integration [START_REF] Lemieux | Monte carlo and quasi-monte carlo sampling[END_REF]. This method has barely been studied for the numerical resolution of thermal radiation. First works on QMC [START_REF] Sarkar | A comparative study of pseudo and quasi random sequences for the solution of integral equations[END_REF][START_REF] Kersch | Radiative heat transfer with quasi-monte carlo methods[END_REF] in radiative heat transfer have highlighted higher accuracy of QMC compared to MC in one-dimensional slabs simulations as well as in a bidimensional problem with radiative heat transfer between surfaces in semi-conductor applications. With one notable exception [START_REF] Marston | Geometric optimization of radiant enclosures containing specularly-reflecting surfaces through quasi-monte carlo simulation[END_REF], the few earlier studies show that this method has not been considered for practical and costly simulations, with modern standards. To the best of our knowledge, no applications of Quasi-Monte Carlo to 3D participating media have been reported in the literature. Additionally, in order to get an error estimation which is, if not necessary, appreciable, a randomization of Quasi-Monte Carlo is needed (Randomized-QMC or RQMC). Only the recent study by Marston et al. [START_REF] Marston | Geometric optimization of radiant enclosures containing specularly-reflecting surfaces through quasi-monte carlo simulation[END_REF] has considered such an RQMC technique in radiative heat transfer simulation: radiative exchanges between surfaces in a transparent medium are solved with a Halton low-discrepancy sequence and a random shift randomization technique. A Randomized-QMC method will also be retained here and the resulting modified statistical estimator for the computed radiative power will be thoroughly characterized. The error estimation provided by RMQC, which is not available with Quasi-Monte Carlo only, also allows a clear and consistent comparison with the classical Monte Carlo method.

The present study focuses on the assessment of randomized Quasi-Monte Carlo methods to solve radiative heat transfer in comparison with MC methods. This is done in several configurations (turbulent channel flow, jet flame, confined premixed flame) with instantaneous 3D solution fields obtained from DNS or LES. The considered cases are characterized by several degrees of complexity: high pressure, presence of soot and finally participating walls. This enables to study the impact of QMC methods in practical 3D cases and to anticipate its benefits in coupled high-fidelity simulations based on LES or DNS.

The principle and properties of Quasi-Monte Carlo method and its randomized variant are illustrated in Sec. 2 on a simple numerical integration problem. Section 3 details the Monte Carlo solver in charge of describing the radiative transfer and computing the radiative power fields with MC or QMC methods. The three application cases are then presented in Sec. 4. The accuracy and efficiency of the retained randomized QMC method are finally studied in Sec. 5.

Monte Carlo and Quasi-Monte Carlo integration

Definition of a simple test case and Monte Carlo integration

In a general multivariate integration problem, the goal of MC methods is to estimate:

I = V f (x)dx ( 1
)
where f is a real-valued function defined over the given state-space volume V . In order to illustrate the properties of MC and QMC methods, the following simple 2-D integral is considered in this section:

I = 2π 0 π 0 F (θ, φ) sin(θ) 2 dθ 1 2π dφ. (2)
I is then expressed as the expectation of F (θ, φ), I = E[F ], based on the probability density functions p θ (θ) = sin(θ)/2 on the interval [ 0 , π ] and the uniform distribution p φ (φ) = 1/(2π) on the interval [ 0 , 2π ]. The Monte Carlo estimate for the integral is then defined as

Q(N ) = 1 N N i=1 F (x i ), (3)
where x i are samples from the random vector x = (θ, φ) whose components follow the probability density functions p θ (θ) and p φ (φ), respectively. The samples x i = (θ i , φ i ) are determined through two independent random numbers R φ and R θ uniform between 0 and 1 as

φ = 2πR φ (4) θ = cos -1 (1 -2R θ ) (5)
The function F (θ, φ) to integrate is chosen as:

F (θ, φ) = sin 2 (φ/2) sin(θ) (6)
and the exact result for the integral in Eq. ( 2) corresponds to I = π/8. The intrinsic standard deviation is defined as σ 2 int = E (F -I) 2 . Its exact value corresponds here to

σ 2 int = 1 4 -π 8 2 ≈ 0.31.
The MC results for the integral defined in Eq. ( 2) and ( 6) are shown for 900 000 realizations in Fig. 1. The convergence Q(N ) towards the expected value π/8 (red dashed line) illustrates the law of large numbers where MC methods are rooted. Additionally, the Central Limit Theorem allows for deriving error estimates of the computed value from the samples standard deviation σ computed as

σ = 1 √ N -1 N i=1 (F (x i ) -Q(N )) 2 1/2 . ( 7
)
The standard deviation σ is an approximate of the intrinsic standard deviation σ int and the probability density function of Q(N ) tends toward a normal distribution centered in I with a standard deviation equal to σ int /N 1/2 ≈ σ/N 1/2 . This property provides straightforwardly confidence intervals for the estimate Q(N ). These confidence intervals are also shown in Fig. 1. They are seen to decrease with N , which allows to control the number of iterations necessary to reach a given accuracy. The error ε(N ) = |Q(N ) -I| is plotted in Fig. 2. Although the curve is polluted by noise inherent to the statistical estimation of the result, the error is seen to decrease with N approximatively at a constant rate in this logarithmic plot. This is expected and can be consistently characterized by considering a statistical moment of the error to get rid of statistical noise. In particular, the standard )] is computed from 100 independent Monte Carlo simulations and is also shown in Fig. 2. After several dozens of realizations, the asymptotic Central Limit Theorem (CLT) becomes valid and σ[Q(N )] closely follows the plotted law σ int N -1/2 . The statistical error of the MC estimate Q(N ) decreases then according to 1/ √ N . Therefore, in order to decrease the error by a factor ten, the number of samples must be increased by a factor 100. This well-known feature of slow convergence is the main drawback of MC methods. Nevertheless, the error estimate,

deviation of the MC estimate σ[Q(N )] = E (Q(N ) -I) 2 1/2 is directly related to the quadratic mean of the error, E ε(N ) 2 1/2 . σ[Q(N
σ[Q(N )] ≈ σ Q,1 (N ) = σN -1/2 , (8)
allows for a control of the accuracy of the computed value I from the sample values, which is very appreciable. Another methodology to estimate the error, equivalent to this first one, consists in breaking up the samples resulting in Q(N ) into M subsets. For each subset i ∈ [1, M ], a MC estimate Q i (P ) can be computed where P = N/M is the number of samples in the subset. The MC estimations from the total number N of samples is simply related to the MC estimates in the subsets as

Q(N ) = 1 M M i=1 Q i (P ). (9)
The M subsamples Q i (P ) are independent estimations of the same quantity Q(P ) whose expectation is I ≈ Q(N ). The standard deviation of estimates Q i (P ) can then be estimated as

σ[Q(P )] 2 ≈ 1 M -1 M i=1 [Q i (P ) -Q(N )] 2 . ( 10
)
The CLT theorem states that σ[Q(N )] 2 ≈ 1 M σ[Q(P )] 2 which finally yields an alternative estimation of the standard deviation of the MC estimate from all samples [START_REF] Modest | Radiative heat transfer[END_REF][START_REF] Lemieux | Monte carlo and quasi-monte carlo sampling[END_REF]:

σ[Q(N )] ≈ σ Q,2 (N ) = 1 M (M -1) M i=1 [Q i (P ) -Q(N )] 2 1/2 (11)
The two estimates for the MC accuracy that are σ Q,1 (N ) and σ Q,2 (N ) are plotted in Fig. 3 as a function of the total number of samples N . For the second estimate, σ Q,2 (N ) is evaluated from a fixed number M = 1 000 of subsets. Results for σ Q,2 (N ) are then available for N = 1 000 (P = 1), N = 2 000 (P = 2), ..., up to N = 1 000 000 corresponding to a number P = 1 00 of samples per subset.

For each N , the MC result Q(N ) is identical. Only the error estimate of σ[Q(N )] differs. Figure 3 shows that both formulations in Eqs. [START_REF] Jones | Combination of dom with les in a gas turbine combustor[END_REF] or [START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF] are equivalent in the retained conditions (large N and large M ). As expected, both MC error estimates decrease with the square root of the number of realizations N , following the law of large numbers. Figure 3 also shows the evolution of the exact MC root-mean-square error (rms) σ[Q(N )], which can be here computed since the reference solution, I, is known. However, in practical applications as the ones presented later, no analytical solution exists. As a consequence, the MC error cannot be computed and the best practice is to compute the estimations of the MC error, σ Q,1 (N ) or σ Q,2 (N ), which become undistinguishable from the error rms for large N as seen in Fig. 3. Such estimates of the accuracy are commonly considered to assess the performance of MC methods.

Quasi-Monte Carlo methods

The accuracy of a MC calculation is proportional to the variance σ 2 int and inversely proportional to the number of samples as N -1/2 . A possibility to accelerate the Monte Carlo convergence rate as a function of N consists in the use of alternative samplings such as low-discrepancy sequences. The resulting quadrature method is called Quasi-Monte Carlo (QMC).

Low-discrepancy sequences

This type of sequences aims at a more uniform filling of the sampled state space. This is illustrated in Fig. 4 where a MC sequence from a pseudo-random number generator is compared to a low-discrepancy Sobol sequence in two dimensions. The clustering seen with classical random sampling leaves regions with few samples, thus making the MC convergence slower compared to the QMC method where points are more uniformly distributed. As the Monte Carlo method, QMC methods are a numerical method to compute high-dimensional integrals. They rely on the Koksma-Hlawka inequality [START_REF] Hlawka | Funktionen von beschränkter variatiou in der theorie der gleichverteilung[END_REF] which states that for any sequence {u n } and any bounded function f defined in [0, 1] d , the integration error is bounded as:

[f ] = |Q(N ) -I(f )| ≤ V [f ]D * n (u 1 , ..., u n ) ( 12
)
where D * n is the discrepancy of the set of points u i 1 ≤ i ≤ n and V [f ] is the variation of f in the sense of Hardy-Krause. The discrepancy D * n measures how much the sequence deviates from uniform filling. Building low-discrepancy sequences can allow then for more accurate cubature methods than relying on random samples of independent points. In QMC methods, random sequences used in MC are replaced by a deterministic sequence whose points are not independent anymore. The Central Limit Theorem doest not apply then to deterministic Quasi-Monte Carlo simulations.

Details about the construction of low-discrepancy sequences are given in [START_REF] Lemieux | Monte carlo and quasi-monte carlo sampling[END_REF]. Several algorithms have been proposed for the generation of such sequences, and the Halton [START_REF] Halton | On the efficiency of certain quasi-random sequences of points in evaluating multidimensional integrals[END_REF], Sobol [START_REF] Sobo | Uniformly distributed sequences with additional uniformity properties[END_REF] and Faure [START_REF] Faure | Discrepancy of sequences associated with a numeration system (in s-dimension)[END_REF] sequences are the most usually considered. Morokoff et al. [START_REF] Morokoff | Quasi-monte carlo integration[END_REF] demonstrated that the Halton sequence generally outperforms the others when the integral dimensionality is lower than six. For high-dimensional integrals, Sobol sequences exhibit better convergence properties than the Faure or

Halton sequences [START_REF] Galanti | Low-discrepancy sequences: Monte carlo simulation of option prices[END_REF]. Since the integration dimension can be very high in the radiative heat transfer involving diffuse wall reflexions, the Sobol sequence has been retained in this study. The implemented construction algorithm follows the algorithm described in Ref. [START_REF] Joe | Constructing sobol sequences with better two-dimensional projections[END_REF].

The concept of discrepancy only provides bounds for the integration error and is in fact of little help in practice: D * n and V (f ) are not always known, and the upper bound actually significantly overestimates the actual error [START_REF] Boyle | Monte carlo methods for security pricing[END_REF][START_REF] Joy | Quasi-monte carlo methods in numerical finance[END_REF]. Besides, low-discrepancy sequences not being based on independent and identically distributed (i.i.d) points, the CLT does not apply and no general theorem on asymptotic convergence rate properties is then available. While being far more accurate than MC methods, relying on a finite-size sequence for QMC integration does not provide any inherent error estimate. A practical method to obtain error estimates for QMC is based on the randomization of the low-discrepancy sequences and the resulting methods are called Randomized Quasi-Monte Carlo (RQMC) [START_REF] L'ecuyer | Recent advances in randomized quasi-monte carlo methods[END_REF].

Randomized Quasi-Monte Carlo methods

The core idea behind RQMC is to randomize the existing low-discrepancy sequences. This can be done by building M low-discrepancy sequences of P points, where each sequence is a randomized version of the original low-discrepancy sequence and is independent from the others. The randomization technique used in this study is discussed later.

In this way, it is possible to create a random sample of M i.i.d. quasi-random estimators Q QMC i (P ) of P points. The RQMC result is built as the average as in Eq. ( 9),

Q RQMC (N ) = 1 M M i=1 Q QMC i (P ) (13)
and its variance is estimated as in [START_REF] Jones | Combination of dom with les in a gas turbine combustor[END_REF] as

σ Q RQMC (N ) 2 ≈ 1 M (M -1) M i=1 Q QMC i (P ) -Q RQMC (N ) 2 . ( 14
)
Contrary to MC methods where two error estimates, σ Q,1 (N ) and σ Q,2 (N ), can be used, only an equivalent to σ Q,2 (N ) can be built for Quasi-Monte Carlo methods.

Many techniques exist to randomize low-discrepancy sequences, and they have to respect two conditions: the first one is that each point in the randomized point set follows a law of uniform probability on the unity hypercube [0, 1] d for the estimator to be not biased, and the second is that the low discrepancy of the new sequence is preserved after the randomization. The second condition guarantees that the fast convergence of the randomized QMC sequence is not penalized.

Shifting, Full Scrambling and Linear Scrambling are such randomization techniques. Randomizations based on a shift (the simplest one is proposed in Ref. [START_REF] Cranley | Randomization of number theoretic methods for multiple integration[END_REF]) do not sufficiently scramble the point set and, for a certain type of functions, they risk to perform worse than MC integration [START_REF] Lemieux | Monte carlo and quasi-monte carlo sampling[END_REF]. This drawback is not present in the full scrambling approach [START_REF] Owen | Scrambled net variance for integrals of smooth functions[END_REF]. However, such a technique is expensive in terms of memory requirement. For this reason a technique of Linear Scrambling, and more precisely the I-binomial scrambling [START_REF] Tezuka | I-binomial scrambling of digital nets and sequences[END_REF], has been retained in this study, since it preserves the convergence results of the Full Scrambling insuring a good compromise between memory requirement and computational cost. The key principle of such scrambling techniques is the application of a permutation to each digit issued from the decomposition of the sequence points in base b (b=2 for the considered Sobol sequence). The I-binomial scrambling is part of the techniques of linear scrambling, where the permutations are constructed through linear combinations. More details on the construction algorithm can be found in Refs. [START_REF] Lemieux | Monte carlo and quasi-monte carlo sampling[END_REF][START_REF] Tezuka | I-binomial scrambling of digital nets and sequences[END_REF].

Application of Quasi-Monte Carlo integration

The characteristics of Quasi-Monte Carlo integration are here illustrated for the computation of the integral in Eq. ( 2) before considering its application to radiative heat transfer in the next sections.

The integral in Eq. ( 2) is first solved with a non-randomized QMC method. In Fig. 5, the evo- expected value is attained for a number of realizations smaller than the one needed by MC, and the absolute error in QMC is smaller than MC for almost every value of N (especially for large values of N). However, low-discrepancy sequences are not based on i.i.d. sample of points, it is then not possible to get error estimates when applying QMC methods. In order to take advantage of both easy error estimation as done in MC and fast convergence of QMC, Randomized QMC methods are applied 220 to the solution of the integral of Eq. ( 2), where random samples of quasi-random estimators are created.

lution of Q(N ) (left)
MC and Randomized QMC methods are now compared consistently by considering a thousand independent integral estimates, i.e. M = 1 000, based on MC and QMC quadrature for sequences of size P . The average of the M samples yield the respective MC and Randomized QMC estimates based on a total of N = M P evaluations. In Fig. 6, the evolution of the quadrature estimate Q(N) is plotted as a function of P for a fixed M . Similarly to standard QMC, the randomized QMC method reaches the exact value with a smaller number of realizations than MC integration. Note that the numerical estimates are already quite accurate for small values of P because of the retained high value of M chosen to momentary get rid of uncertainties due to the selected value. In addition to fast convergence attributed to standard QMC, the randomized Quasi-Monte Carlo method benefits from an error estimate. Contrary to classical Monte Carlo integration which provides two error estimates through Eq. ( 8) or [START_REF] Jones | Combination of dom with les in a gas turbine combustor[END_REF] , RQMC relies on Eq. ( 14) which is only equivalent to Eq. ( 8). The error estimate σ[Q(N )] for RQMC method is compared to the corresponding Monte Carlo result in Fig. 7 where P is varied while M is still fixed to 1 000. As expected from theory, the MC convergence rate (dashed yellow line) scales as N -1/2 ∝ P -1/2 . For the investigated case, the upper bound of the RQMC error (dashed red line) appears to scale as P -1 . The error estimate confirms the significantly enhanced convergence rate of the RQMC method. A standard deviation of 10 -4 is obtained for P ≈ 300 whereas it is not yet achieved with standard MC with P = 10 4 . This estimate of the RQMC error allows for controlling the number of realizations needed to reach a desired accuracy.

Once again, the estimated error is already small for a small P value because of the chosen value M . Looking at the lower plot of Fig. 7, it is possible to observe the presence of several sub-peaks in σ Q RQMC (N ) (blue solide line). These peaks correspond to those realizations where P = {2 2 , 2 3 , ... , 2 13 } and are characterized by a higher accuracy. This is a property of the Sobol sequence which performs better for powers of two. Thus, one can define two laws for the convergence rate that bound the error estimate: an upper bound which was already introduced and a lower bound that correspond to the point P = 2 i . Both are drawn in Fig. 7. For the studied example, the highest convergence rate obtained for powers of 2 scales as P -1.27 . For Quasi-Monte Carlo simulations with a fixed number of realizations, it is then optimal to choose it as a power of 2.

The dependency of RQMC results on the number M of random QMC estimate is now considered with a fixed size of the sequences set to P = 1 000. Figure 8 presents the error estimate σ[Q(N )] as a function of the number of trials, M , for both MC and RQMC. As expected from the CLT which states

that σ[Q(N )] 2 ≈ σ[Q(P )] 2 /M
, both MC and RQMC errors scale as M -1/2 . The shift between both curves is due to the higher accuracy of RQMC compared to MC for a given sequence size P .

The computational cost of each method depends on the total number of realizations N , i.e. the product of both M and P . Previous results have shown that MC and RQMC convergence rates scale as M -1/2 P -β where β = 1 2 in the MC method, while it is appreciably larger in RQMC and unknown a priori for a given case. In standard MC simulation, the convergence rate depends then eventually as N -1/2 and is not sensitive to the M -P distribution. The number of MC trials M used to compute the error estimate must nonetheless be large enough for Eq. ( 11) to be reliable. A minimum value of M often admitted in the literature is 10 [START_REF] Lemieux | Monte carlo and quasi-monte carlo sampling[END_REF]. In RQMC, the trade-off between M and P for a given

N is not at all trivial. Error estimates as a function of the total number of samples N are shown in Fig. 9 for both MC and RQMC computations with different M -P repartitions. On the one hand, all MC results are similar and tend to the same asymptotic convergence rate as previously explained.

On the other hand, it can be clearly seen that each M -P combination yields a different convergence rate for RQMC. For a same computational cost (identical N ), the combination with small values of M achieves significant higher accuracy. This is due to the M -1/2 P -β convergence law with β > 1/2: it is more interesting to increase the size of the low-discrepancy sequence, P , for a same N in order to benefit from the QMC enhancement. The choice of M is then critical in RQMC method.

However, M cannot be too small, otherwise the error estimate in Eq. ( 14) can not be trusted. A compromise must then be found. A sensitivity analysis has been performed to study the evolution of previously. This value of M is retained for the numerical setup of the following simulations.

P )] = 1 M -1 M i=1 Q QMC i (P ) -Q RQMC (N ) 2 1/2
as a function of the number of trials M and for P = 10 000.

Monte Carlo numerical solver of radiative heat transfer

The in-house code Rainier [START_REF] Zhang | Optimized emission-based reciprocity monte carlo method to speed up computation in complex systems[END_REF][START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF][START_REF] Koren | Multiphysics simulation combining large-eddy simulation, wall heat conduction and radiative energy transfer to predict wall temperature induced by a confined premixed swirling flame[END_REF][START_REF] Palluotto | Comparison of monte carlo methods efficiency to solve radiative energy transfer in high fidelity unsteady 3d simulations[END_REF] is described in the present section. It solves the radiative transfer equation (RTE) with a Monte Carlo approach while relying on standard Monte Carlo or Randomized Quasi-Monte Carlo integration. Details of the solver, its adaptation to QMC integration and the retained description of radiative properties are given in the following.

Ray tracing in the Monte Carlo solver and adaptation to QMC

In the Rainier solver, the local radiative power is computed at each vertex of the computational domain (those where the solution is desired). Then, the departing position of a photon bundle is fixed and it is not determined by a random variable. Thus, a ray tracing starting at a given point is initially characterized by a direction determined the unit vector denoted as ∆, or equivalently by two angles θ and φ, and a frequency ν generated according to a given probability density function depending on the resolution method.

A ray is tracked along its optical path with possible absorption by the participating medium or walls, wall reflexions and exiting of the computational domain. To keep the simulation affordable, the ray tracing is interrupted when the equivalent transmission from the departing point is below a predefined threshold τ min . In other words, a ray is stopped when its remaining energy has been lowered by a factor corresponding to (1τ min ). Here, τ min is set equal to 0.01 in all the investigated cases.

During its optical path, a ray hitting a wall may be reflected. For black walls or specular reflexion, no additional random number than the previous three for frequency and departing direction is necessary.

The corresponding three-dimensional space is then sampled with either a pseudo-random generator (the one used in the present work is L'Ecuyer's Multiple Recursive Generator MRG32k3a [START_REF]Good parameters and implementations for combined multiple recursive random number generators[END_REF]) or a three-dimensional Sobol sequence for QMC integration. However, cases with diffuse reflexions as later considered in confined configurations require a specific attention. In standard MC, the reflected direction then needs to be randomly generated and the integration dimension d of the problem increases as d = 3 + 2r where r is the number of reflections performed. The number of reflexions that one ray can undergo is then unknown, making d undetermined as well. While this is transparent for standard MC where the pseudo-random generator provided i.i.d samples, this is critical for QMC integration that requires to a priori know the dimension of the integration problem to build the lowdiscrepancy sequence. In the present integration of such sequences, the maximum number of reflexions r max is determined from the worst case scenario and is used to build a d-dimensional sequence with

d = 3 + 2r max .
The worst case is considered as a transparent medium with wall diffuse reflectivities equal to 1min , where min is the the lowest wall emissivity considered in the problem setup. After r reflexions, the energy of the ray is scaled down by the factor (1min ) r . Given the criterion τ min introduced earlier to stop the ray, the maximum number of reflexions follows

(1 -min ) rmax < τ min (15)
Thus,

r max > ln(τ min ) ln(1 -min ) (16)
and the dimension d becomes:

d = 3 + 2 ln(τ min ) ln(1 -min ) . ( 17
)
The formula obviously falls apart for purely reflective surfaces ( min = 0), which is unrealistic in practical applications of radiative heat transfer.

Finally, the number of rays issued from a given point of interest is controlled by the prescribed accuracy. The standard Monte Carlo accuracy is estimated by assembling rays in several batches (M = 20) and using the estimate in Eq. ( 11) for the quantity of interest (radiative power of flux). This makes the implementation of the Randomized-Quasi Monte Carlo error estimate given by Eq. ( 14)

straightforward. For both methods, the control is carried out on the relative and the absolute value of the standard deviation of the results estimate. Hence, the relative standard deviation of the radiative power is defined as the ratio of the local standard deviation to the local radiative power. However, this relative standard deviation is not appropriate in regions with negligible radiative contribution.

A second control is then considered on the absolute value of the standard deviation which is checked to be lower than a prescribed maximum, typically set equal to a given percentage of the maximum radiative power. When error control is enabled, the prescribed tolerances for the absolute and relative standard deviation allow for a local adaptation of the number of rays et each point.

Emission-based Reciprocity Method

Two backward Monte Carlo methods are implemented in the solver: the Emission Based Reciprocity method (ERM) [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF] and the Optimized-ERM (OERM) [START_REF] Zhang | Optimized emission-based reciprocity monte carlo method to speed up computation in complex systems[END_REF]. The main interest in such methods is the possibility to dynamically and independently control the local convergence of the computed radiative fields. Photon bundles are traced from volume or surface points of the computational domain where the radiative power or wall flux is predicted. The independent treatment of points of interest allows for a high degree of scalability. A scalability test has been performed on a cluster equipped with Intel E5-2690 processors on a configuration characterized by a computational domain of 8 millions cells.

The number of cpu cores is varied from 120 up to 1920. Results, displayed in Fig. 11, show that a very good scalability is obtained. Besides, ERM and OERM are reciprocal MC methods that exactly enforce the Reciprocity Principle [START_REF] Case | Transfer problems and the reciprocity principle[END_REF], that is only statistically fulfilled by other Monte Carlo methods. In a reciprocal formulation the radiative power per unit volume of the node i of the discretization domain is calculated as the sum of the exchanged powers P exch ij between the node i and all the other cells j crossed or encountered by the optical paths generated by i.

P i = N j=1 P exch ij ( 18
)
where N is the number of volume cells and faces in which the domain is discretized.

P exch ij
is the difference between the radiative power emitted by j, transmitted by the medium and absorbed by i, and the radiative power emitted by i, transmitted by the medium and absorbed by j:

P exch ij = ∞ 0 κ ν (T i )[I 0 ν (T j ) -I 0 ν (T i )] 4π A ijν dΩ i dν (19)
where κ ν is the spectral absorption coefficient, I 0 ν the equilibrium spectral intensity, Ω i is the solid angle issued from i, ν the photon frequency and A ijν accounts for all the paths between emission from i and absorption in j after transmission and possible walls reflections along the paths. Its expression is detailed in [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF].

Introducing the volumetric power P e i (T i ):

P e i (T i ) = 4π ∞ 0 κ ν (T i )I 0 ν (T i )dν (20)
Equation ( 19) can be written as:

P exch ij = P e i (T i ) ∞ 0 I 0 ν (T j ) I 0 ν (T i ) -1 4π A ijν f i (∆, ν)dΩ i dν (21)
where f i (∆, ν) is the joint PDF with ∆ the direction.

f i (∆, ν) can be separated in two independent parts:

f i (∆, ν)dΩ i dν = f ∆i (∆)dΩ i f νi (ν)dν (22)
with:

f ∆i (∆) = 1 4π ; f νi (ν) = κ ν (T i )I 0 ν (T i ) ∞ 0 κ ν (T i )I 0 ν (T i ) (23)
f νi (ν) is the frequency distribution function and it is equal to the emission distribution function at the temperature T i of the emitting node. Then, in the ERM formalism, the frequency ν n associated with a realization n is obtained by a uniform random number R n between 0 and 1, and ν n is determined from the equation:

R n = νn 0 f νi (ν)dν = νn 0 κ ν (T i )I 0 ν (T i )dν ∞ 0 κ ν (T i )I 0 ν (T i )dν (24)
In the OERM method, the frequency distribution function f νi is based on the emission distribution at the maximum temperature encountered in the system, T max ; as a consequence, ν n in the OERM formalism is determined from:

R n = νn 0 κ ν (T max )I 0 ν (T max )dν ∞ 0 κ ν (T max )I 0 ν (T max )dν (25)
The direction ∆ is determined in spherical coordinates by the polar angle θ and the azimuthal angle φ, which are determined according to Eqs. [START_REF] Habibi | Impact of radiation models in cfd simulations of steam cracking furnaces[END_REF][START_REF] Coelho | Spectral radiative effects and turbulence/radiation interaction in a non-luminous turbulent jet diffusion flame[END_REF]. Statistical estimation of the radiative power P i at the node i, indicated as Pi , can be obtained by summing all the contributions of the N i optical shots traced from i. And its final expression, as it is computed in the presented simulations, is:

Pi = P e i N i Ni n=1 Mn m=1 I 0 νn (T m ) I 0 νn (T i ) -1 τ νn α n,m,νn (26)
where m = 1 is the cell i, while m = M n represents the last cell crossed by the n th optical path originating from i. Then, α n,m,νn is the spectral absorptivity of the n th optical path in the m th cell crossing, while τ νn is the spectral transmissivity from i to the cell m accounting for eventual wall reflections.

From the Eq. ( 26) it is possible to see that the radiative power exchanged between two cells at identical temperatures is rigorously null, leading to more accurate results compared to conventional formalisms.

Gas radiative properties

Concerning the gas radiative properties, a good compromise between accuracy and computational cost is found with a narrow-band approach: the correlatedk or ck model by [START_REF] Goody | Atmospheric radiation: theoretical basis[END_REF] based on updated parameters due to [START_REF] Rivière | Updated band model parameters for h2o, co2, ch4 and co radiation at high temperature[END_REF], which have been obtained for the range of temperature and pressure of interest, by using the CDSD-4000 database for the absorption spectra of CO 2 [START_REF] Tashkun | Cdsd-4000: High-resolution, high-temperature carbon dioxide spectroscopic databank[END_REF]

Validation of the RQMC implementation

In order to validate the RQMC implemented in the solver Rainier, a 1D slab geometry, for which a semi-analytical solution of the radiative transfer equation can be obtained from exponential integral functions [START_REF] Modest | Radiative heat transfer[END_REF], is retained and two test cases are performed. The radiative power per unit volume is computed within a slab normal to X and bounded by two gray and opaque infinite planes. The slab thickness is L = 0.1 m and the walls are characterized by a diffuse emissivity of 0.6 and a wall temperature T w = 500K. In the first test, denoted as Case a, a gray gas (constant κ ν ) with an optical thickness equal to 1 is considered. The second test, called Case b, represents more realistic conditions: a medium with an homogeneous composition (CO 2 and H 2 O molar fractions are, respectively, 0.116 and 0.155) is considered. Spectral radiative properties of the gaseous mixture are taken into account through the aforementioned c-k model. For both investigated cases, the gas is at atmospheric pressure and has a parabolic temperature profile T (X) = aX 2 + b, where X is the normal to the walls, b is the maximum temperature, equal to 2500 K, and a = Tw-b (L/2) 2 . For each investigated case, MC and RQMC simulations are performed by imposing convergence criteria, which allows to control the accuracy of the solution and have comparable profiles between MC and RQMC fields.. Profiles of radiative power across the normal to the slab, X are presented in Fig. 12, for Case a and in Fig. 13 for Case b. Black solid lines represent the semi-analytical reference solution in both cases. Figures 12 and13 show that, for both investigated cases, the MC and RQMC results overlap on the reference profile. This validates the implementation of RQMC in the Rainier solver. In the next section, more practical configurations are retained for the thorough comparison of MC and RQMC methods.

Investigated configurations

Three 3D configurations are investigated: a canonical flow configuration (turbulent channel flow at high pressure) and two combustion applications (a turbulent sooting jet flame and a combustion chamber featuring a premixed swirled flame). These three cases are characterized by different degrees of complexity and enable to study the impact of the radiative heat transfer in different practical situations: at high pressure, in reacting flows, in presence of soot or walls. This allows for comparing MC and RQMC methods behavior in several contexts. Besides, all temperature and composition fields are extracted from direct numerical simulations and large-eddy simulations. The numerical methods 410 are then assessed in configurations relevant to coupled high-fidelity simulations.

Case 1: Turbulent channel flow

The first investigated application is a turbulent channel flow (case C3R1 from [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]) whose dimensions are: L x × L y × L z = 0.62 m × 0.2 m × 0.3 m. The channel is filled by an homogeneous gas mixture of CO 2 -H 2 O -N 2 with a corresponding molar fraction of 0.116 -0.155 -0.729. Its pressure is 40 bars and two iso-thermal walls (at 2050 K on the top and 950 K on the bottom) with a total hemispherical emissivity of 0.8 confine the flow. Its structured computational grid is made by 4.2 millions nodes.

In Fig. 14, snapshots of 3D direct numerical simulations taken from [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF] are shown: on the top the instantaneous temperature field and on the bottom the radiative power field are given. At the given high pressure, the optical thickness of the medium is large. The effect of a large optical thickness can be seen in the very thin layer close to the wall: the bottom one is dominated by the absorption, while the one close to the hot wall is dominated by emission. Consequently there is a negligible radiative interaction between the two walls, and the whole domain is characterized by absorbing (yellow pockets) and emitting (blue pockets) regions. The emitted thermal radiation is absorbed in the vicinity of the emitting point.

Case 2: Non-confined sooting jet flame

The second configuration is the turbulent jet diffusion flame experimentally studied at Sandia [START_REF] Zhang | Design of "model-friendly" turbulent non-premixed jet burners for c2+ hydrocarbon fuels[END_REF].

It is fed with pure ethylene at 300 K and is stabilized with the help of pilot flames fed by ethylene and air. Unlike the previous configuration, the participating medium is heterogeneous and contains soot particles. The presence of soot particles make the radiation heat transfer very important.

In Fig. 15, soot volume fraction and temperature fields extracted from 3D Large Eddy Simulations of [START_REF] Rodrigues | Coupling an les approach and a soot sectional model for the study of sooting turbulent non-premixed flames[END_REF] are shown along with the computed instantaneous radiative power field (left). Soot particles are present in richer regions of the flame and maximum instantaneous soot volume fraction is nearly 10 ppm. Radiative power issued from soot particles is located also in richer regions dominated by emission. Similarly, the radiative power is also negative in hot near stoichiometry burnt gases. On the other hand, lean regions where temperature is lower than 1000 K are dominated by the absorption from CO 2 and H 2 O.

When soot radiation is included in the calculation, their absorption spectra need to be taken into account; then 93 spectral bands are introduced between 150 and 29 000 cm -1 , of which 44 overlap with the gas bands. The soot absorption coefficient is modeled through the Rayleigh's theory [START_REF] Modest | Radiative heat transfer[END_REF]:

κ soot ν = C 0 f v ν with C 0 = 36πnk (n 2 -k 2 +2) 2 +4n 2 k 2
where n and k are the real and imaginary part of the complex index of refraction of soot particles, taken as equal to m = nik = 1.57 -0.56i [START_REF] Smyth | The elusive history of m = 1.57-0.56 i for the refractive index of soot[END_REF], ν is the wavenumber (m -1 ) and f v is the soot volume fraction. Scattering by the soot particles is here neglected.

Case 3: Combustion chamber

The third investigated application is a laboratory scale combustor (30 cm long with a square section of 81 cm 2 ) studied experimentally [START_REF] Guiberti | Analysis of topology transitions of swirl flames interacting with the combustor side wall[END_REF][START_REF] Guiberti | Impact of heat loss and hydrogen enrichment on the shape of confined swirling flames[END_REF] and numerically [START_REF] Mercier | Experimental and numerical investigation of the influence of thermal boundary conditions on premixed swirling flame stabilization[END_REF][START_REF] Koren | Multiphysics simulation combining large-eddy simulation, wall heat conduction and radiative energy transfer to predict wall temperature induced by a confined premixed swirling flame[END_REF]. The chamber is operated at atmospheric pressure and a turbulent premixed flame of CH 4 , H 2 and Air is stabilized. The flame is characterized by a thermal power of 4 kW and is confined by cold quartz windows. Such configuration is typical of industrial furnaces. Figure 16 shows the temperature field extracted from Large Eddy Simulations [START_REF] Koren | High-fidelity multiphysics simulation of a confined premixed swirling flame combining large-eddy simulation, wall heat conduction and radiative energy transfer[END_REF] on a longitudinal section of the chamber along with the presently computed radiative power field. The 3D computational grid is made of 1.26 millions points.

As in the previous configuration, heterogeneities of temperature and of participating species are present. The complexities seen in the previous cases (heterogeneities and participating walls) are simultaneously taken into account in this configuration and they represent typical situations encountered in a real enclosed application. The temperature field (Fig. 16 top) shows that burnt gases issued from the flame are hotter than the ones close to the walls. Because of the presence of recirculation zones between the flame and the walls, hot gases are recirculated toward the chamber bottom while they are cooled down by wall convective heat transfer. The radiative power field indicates that cold gases close to the walls are dominated by absorption (red regions/positive radiative power), while the higher temperature regions mostly emit radiation (negative values of radiative power). In the present radiative heat transfer computations, the walls are considered as opaque with a total hemispherical emissivity of 0.75. imposing an identical number of rays at all the computed points, while no convergence criterion is imposed. Such an analysis allows to evaluate the level of convergence, i.e. the standard deviation achieved at every node of the computational domain where the solution is desired. On the contrary, in the second test, Test 2, two convergence criteria are imposed: relative and absolute standard deviation.

In this case, the number of rays is not set a-priori but it spatially varies over the computational domain.

Table 1 summarizes the setup used for the tests performed on each application.

For all the computations, inlet and outlet are considered as non-reflecting ( = 1) with a far-field temperature of 300 K (except for Case 1, where periodic conditions are considered). The retained method is the ERM method [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF] in its MC and RQMC version.

Results: three practical applications

For each configuration, results of radiative heat transfer simulations are shown. Tests with a fixed number of rays, also referred as realizations, are first carried out before considering tests with controlled convergence.

Tests with a fixed number of realizations

MC and RQMC methods are compared in simulations where the number of realizations is imposed (Test 1 in Tab. 1) in order to compare the local convergence attained in all the points of the computational grid. The accuracy of the computations is quantified as the standard deviation of the estimated radiative power following Eqs. 11 and 14. In Fig. 17, the absolute standard deviation scaled by the the maximum radiative power in the domain is shown for Case 1 for MC (top) and RQMC simulations (bottom). With the same number of realizations, the RQMC computation allows for converging all the points of the computational domain with a significantly smaller absolute standard deviation compared to the classical MC approach. For both methods, the resulting standard deviation field is not uniform.

In particular, the standard deviation is larger close to the cold wall. This behavior is expected with the ERM approach whose drawback is a poorer convergence in cold regions with significant absorption. Additional tests are performed in order to check the convergence rate of MC and RQMC simulations as a function of the number of realizations N for the three investigated cases (M is still fixed to 20).

For each configuration, the relative standard deviation is measured in several characteristic points of the domain, highlighted in Figs. dashed lines) and RQMC (gray dotted line) are also highlighted on the figures. In the case of RQMC computations where no theoretical results provide the asymptotic convergence rate, the corresponding 510 lines illustrate the convergence rate of the method by considering a slope computed by a linear fitting of the RQMC standard deviation on the different probes. The proportionality factor is however different and all MC results appear as translated from each other in the logarithmic plots. As indicated by Eq. 8, this is the result of different intrinsic standard deviation (estimated as σ) at the considered points. σ depends on the specific radiative exchanges for the given location and on the Monte-Carlo algorithm solving for radiative heat transfer. The ERM methods is used here, which is known to converge more poorly in colder regions as confirmed in Figs. 20 and21. The figures reveal that, while RQMC simulations are more accurate than MC ones and present slightly different slopes, the RQMC convergence laws also exhibit translation shifts between the different points. Although the observed distance between RQMC results is different from the one seen with MC results, the same ranking of accuracy between the different probes is kept for both approaches. The RQMC results are then sensitive to the intrinsic standard deviation of the integral to compute.

Figure 22 shows the standard deviation achieved in MC and RQMC simulations of the combustion chamber (Case 3), for the three considered points, F , G and H. Concerning the cold region of burnt gases represented by probe G, the RQMC standard deviation is higher than MC one for N < 10 3 , showing again that RQMC does not always guarantee more accurate results with a reduced number of evaluations. Nonetheless, for higher values of N , the asymptotic stronger convergence law of RQMC at point G eventually yields the best results. For points F and H at high and intermediate temperature respectively, the RQMC method behaves better than MC for every N . The dotted lines, obtained by fitting the points from N = 10 4 to N = 10 6 , show that RQMC convergence rate scales as N -0.725 for Point F , as N -0.632 for G and as N -0.813 for H. Also in this configuration, the magnitude of the RQMC exponent b remains higher than b = 1/2 for MC methods for the three investigated regions of the combustion chamber.

Controlled convergence tests

In Figs. 20, 21 and 22, the smallest RQMC standard deviation is obtained in points characterized by the highest temperature, A, D and F . Instead, the largest standard deviation is found in the Point B at intermediate temperature for Case 1 and in the colder regions of the domain, Points E and G, for Cases 2 and 3. Considering a uniform number of rays for each point is then not ideal as some locations can be regarded as too accurate while the standard deviation at others points can remain unsatisfactory and suffer from the retained number of realizations N . In practice, a variable field N (x, y, z) is preferred to spatially control the accuracy of the result and yield a more efficient computation of the radiative fields on the whole domain for a given computational cost. MC and RQMC results are then compared in such a situation described by Test 2 in Tab. 1.

Figure 23 shows the field of the number of rays that are tracked independently at each point to perform controlled-convergence simulation using a MC (top) or RQMC (bottom) method for the channel flow application (Case 1). In MC simulations (on the top), it is possible to notice the presence of a region characterized by a large number of realizations that corresponds to the colder region of the channel. This is consistent with the previous observation in Fig. 20: convergence is more difficult to achieve in the lower part of the channel compared to the hotter parts of the configuration. As expected from the enhanced convergence rate of the RQMC method, simulations performed with RQMC (bottom) needs fewer rays compared to MC in order to achieve the fixed accuracy criteria. This is seen in all regions including where the convergence is harder at the bottom of the channel.

The reduced number of tracked rays leads to a lower computational time for the simulation, which is determined in the next section.

The same test is performed for Case 2 where, unlike Case 1, the composition is heterogeneous and the effect of the soot particles on the radiative heat transfer is considered. MC and RQMC methods are compared in Fig. 24. In both cases, regions needing more realizations are the ones characterized by the presence of cooler pockets of burnt gases without soot. It clearly emerges that RQMC outperforms MC also in the second investigated configuration.

Case 3 represents a confined laboratory scale burner and accounts for heterogeneities of species and temperature fields but also for the presence of participating walls. Fig. 25 compares MC and RQMC methods: regions where a large number of rays is needed in order to achieve the convergence are the ones characterized by a strong gradient of temperature in the region between the burnt gases generated by the flame (dominated by emission) and these gases later trapped in the cold outer recirculation zone (dominated by absorption). The difficulty is associated to the small magnitude of the radiative power surrounded by positive and negative regions, making the fulfillment of the criterion on relative standard deviation hard. Nonetheless, the number of realizations is drastically reduced with RQMC in all regions.

CPU Efficiency of Monte Carlo and Randomized quasi-Monte Carlo methods

A more complete comparison can be done by evaluating the efficiency of both MC and RQMC methods. The MC efficiency is evaluated as a function of variance and computational time [START_REF] Lemieux | Monte carlo and quasi-monte carlo sampling[END_REF]. In the present analysis, the local efficiency η i is computed from the variance associated to the point i, . This last term is computed as the number of intersections that a ray traced from point i experiences along its optical path, nb int,i , multiplied by the cost of one intersection, T CP U /nb int,tot :

η i = 1 σ 2 i [Q(N )] T CP U i (27)
In Fig. 26, the ratio of the local efficiencies of RQMC and MC algorithms is shown on longitudinal planes of the three retained applications. For each configuration, almost the whole domain is characterized by an efficiency ratio larger than unity, meaning that the RQMC methods improve the computational efficiency by a factor that spatially varies in the domain and that may get values higher than 3.

In Case 1 (top), characterized by high pressure and homogeneous mixture, it can be observed that the improvement is significant close to the cold wall of the channel, where the convergence is harder to achieve, and on the region at intermediate temperature in the middle of the channel. In Cases 2 and 3, the improvement is obtained almost everywhere in the computational domain.

It is straightforward to deduce that, for a fixed variance, the larger the efficiency is, the smaller the computational time will be. application, 1.8 for the turbulent sooting jet flame and 2.5 for the combustion chamber. Finally, let us note that the obtained speed-up factor associated to convergence laws approximatively scaling as P -1 is similar to the results from [START_REF] Marston | Geometric optimization of radiant enclosures containing specularly-reflecting surfaces through quasi-monte carlo simulation[END_REF] involving surface radiative transfer.

RQMC combined with importance sampling

Previous results have shown that classical MC results can converge faster at different points depending on the corresponding intrinsic variance. Several strategies for variance reduction, such as importance sampling, are common to accelerate the convergence of Monte Carlo simulations. The RQMC results being sensitive to variance effects, it can also benefit from the same strategies derived for MC computations. In fact, the pseudo-random number generator or retained low-discrepancy sequence are set independently of any variance reduction technique. Quasi-Monte Carlo cubature and variance reduction can then be combined to achieve an even more efficient computation. This is here illustrated with the OERM method [START_REF] Zhang | Optimized emission-based reciprocity monte carlo method to speed up computation in complex systems[END_REF] which is an importance sampling technique to improve the poor convergence of ERM in cold regions dominated by absorption.

The OERM method, shortly described in Sec. convergence rate of RQMC-OERM.

Conclusion

Nowadays, coupled simulations involving combustion and radiative heat transfer are more and more used and targeted. Thanks to the increase in computing power and advances in numerical algorithms and solver scalability, the use of accurate Monte Carlo methods in 3D unsteady simulations, such as Direct Numerical Simulations and Large Eddy Simulations, has become affordable. However, such coupled and multi-physics simulations remain very costly, and additional efforts are necessary to make MC methods more efficient. In this study, a technique to improve MC methods is investigated: Quasi-Monte Carlo method, where pseudo-random sequences of MC are replaced by quasi-random sequences also known as low-discrepancy sequences. Sobol sequences are retained in this study. Their advantage lies in a higher convergence rate compared to MC methods, on the other hand their deterministic feature does not make possible to get an error estimation. Thus, a randomization of QMC (RQMC) is needed and among the existing randomization techniques, the one retained in this study is the I-binomial scrambling. The features of MC and RQMC methods are first shown on a 2-dimensional integration problem, highlighting the higher convergence rate of RQMC when compared to MC. The implementation of the RQMC method is then combined to the ERM formalism in the radiative heat transfer solver. After being validated on 1D test cases, it is applied successively to several applications of radiative heat transfer problems. Three configurations are retained, very different from each others:

an academic test case of a channel flow at 40 bar, a non-confined turbulent sooting jet flame and a semi-industrial combustion chamber operating at atmospheric pressure. Simulations results show a significant improvement of RQMC compared to MC method in terms of computational efficiency, eventually leading to a reduction of the computational time of a radiative heat transfer simulation of all the investigated cases by a factor that varies from 1.8 to 2.5, depending on the configuration.

Convergence enhancement is due the stronger decay of the standard deviation with the number of realizations which overcomes the classical MC convergence law. Finally, RQMC is combined to a technique of importance sampling: it is shown that the efficiency of RQMC method can be further improved when it is combined with a variance reduction technique, such as the OERM method considered here for illustration.

Figure 1 :

 1 Figure 1: Evolution of the MC result Q(N) as a function of the number of realizations N . Red horizontal dashed line: the exact value of the integral. Confidence intervals at 90%, 95% and 99% are denoted by yellow, green and blue dashed lines, respectively.

Figure 2 :

 2 Figure 2: Evolution of deterministic error ε(N ) (blue plain line) and mean quadratic error σ[Q(N )] (black plain line) as a function of the number of realizations N . The red dashed line is the theoretical convergence rate σ int N -1/2 .

Figure 3 :

 3 Figure 3: Evolution of the standard deviation as a function of the total number of realizations N in a log-log scale in the two analyzed cases : red solid line is σ Q,1 (N ) from Eq. (8) in 1 MC trial with up to 1 million realizations; yellow dashed line is σ Q,2 (N ) from Eq. (11) with a fixed number M = 1000 of MC trials with up to P = 1000 realizations each; black dashed line is the line of slope -1/2; blue solid line with white circles represent the standard deviation of MC estimate σ[Q(N )].

Figure 4 :

 4 Figure 4: Comparison of the sampling of polar (θ) and azimuthal angle (φ) generated accordingly to the probability density functions defined in Eqs. (4-5) using a a classic pseudo-random sequence (left) and a Sobol sequence (right).

  and of the absolute error ε(N ) = |Q(N ) -I(f )| (right) is shown as a function of the number of realizations for QMC and compared to results obtained with MC estimation. In both plots of Fig. 5, it is pointed out that QMC converges faster than MC: in QMC calculation the

Figure 5 :

 5 Figure 5: Evolution of Q(N ) (top) and error ε(N ) (bottom) as a function of the number of realizations in a log-log plot. Black solid line: standard Monte Carlo (MC); blue solid line: Quasi-Monte Carlo (QMC); red dashed line: expected value.

Figure 6 :

 6 Figure 6: Evolution of Q(N ) as a function of the number of realizations in each trial, P , where N = M P and M = 1 000. Black solid line: Monte Carlo (MC) results. Blue solid line: Randomized Quasi-Monte Carlo (RQMC) results. Dashed red line: exact integration result.

Figure 7 :

 7 Figure 7: Evolution of error estimate computed from standard deviation formula in Eqs. (11) and (14) for MC and RQMC, respectively. Results are reported as a function of the number of realizations in each trial, P , for a fixed number M = 1 000 of integration trials. Black solid line: Monte Carlo results. Blue solid line: Randomized Quasi-Monte Carlo results. Red dashed line: upper bound RQMC convergence rate. Magenta dashed line: convergence rate of RQMC for P = {2 2 , 2 3 , ... , 2 13 }.

Figure 8 :

 8 Figure 8: Evolution of error estimate as a function of the number of trials, M for Monte Carlo (black) and Randomized Quasi-Monte Carlo (blue). Each trial contains P = 1000 points.

Figure 9 :

 9 Figure 9: Evolution of error estimate as a function of the total number of realizations, N = M P , for several combinations of M and P identified by colors. Dashed lines: MC results pointed out by black arrow. Solid lines: RQMC results.

Figure 10 :

 10 Figure 10: Evolution of σ[Q(P )] =

Figure 11 :

 11 Figure 11: Scalability of the Rainier solver with number of cores. Red circles correspond to the performed tests, dashed line corresponds to the ideal performance.

Figure 12 :

 12 Figure 12: Profile of radiative power in Case a across the slab coordinate X computed by MC (blue) and RQMC (red), compared to the reference solution (black solid line).

Figure 13 :

 13 Figure 13: Profile of radiative power in Case b across the slab coordinate X computed by MC (blue) and RQMC (red), compared to the reference solution (black solid line).

Figure 14 :

 14 Figure 14: Instantaneous fields of temperature (top) and radiative power (bottom) on a longitudinal section of the channel. Points A, B and C, whose results are analyzed in Sec. 5, are indicated by black circles on the top figure.

Figure 15 :

 15 Figure 15: Instantaneous fields of radiative power (left), soot volume fraction (center) and temperature (right) on a longitudinal section of the jet flame. Color scales are linear for radiative power and temperature and logarithmic for soot volume fraction. Points D and E, whose results are analyzed in Sec. 5, are indicated by black circles on the right figure.

Figure 16 :

 16 Figure 16: Instantaneous fields of temperature (top) and radiative power (bottom) on a longitudinal section of the burner. Points F , G and H, whose results are analyzed in Sec. 5, are indicated by black circles on the top figure.

4

 4

  The same observation can be done for Cases 2 and 3 shown in Figs.[START_REF] Walters | Rigorous development for radiation heat transfer in nonhomogeneous absorbing, emitting and scattering media[END_REF] and Fig.19in terms of the local relative standard deviation, i.e the ratio between the standard deviation and the local value of the radiative power. It can be seen that RQMC simulations (on the bottom) are more accurate than MC ones (on the top) since, with RQMC, most of the points of the domain achieve a relative standard deviation lower than the one obtained when MC is used.

Figure 17 :

 17 Figure 17: Field of absolute standard deviation scaled by the the maximum radiative power in the domain. The computation is carrie out following Test 1 with MC (top) and RQMC (bottom) in the channel flow configuration (Case 1).

  14, 15 and 16 and identified by letters. In Case 1, three points, A, B and C, are retained: they are located in proximity of the hot wall, at the half-height and close to the cold wall of the channel, respectively. In Case 2, two points, D and E, are chosen to represent, respectively, a hot and rich mixture characterized by a high soot concentration and downstream burnt gases at lower temperature where no soot is present. In Case 3, the probes, F , G and H, are located in the flame region, in the cold outer recirculation zone (ORZ) close to the chamber bottom and further downstream in the ORZ at an intermediate temperature, respectively.

Figures 20 ,

 20 Figures 20,[START_REF] Delatorre | Monte carlo advances and concentrated solar applications[END_REF] and 22 display the relative standard deviation obtained by MC (circles) and RQMC (diamonds) for each studied points in the different configurations. The convergence rates of MC (grey

Figure 18 :

 18 Figure 18: Field of relative standard deviation achieved in Test 1 with MC (top) and RQMC (bottom) in a sooted jet flame application (Case 2).

Figure 19 :

 19 Figure 19: Field of relative standard deviation achieved in Test 1 with MC (top) and RQMC (bottom) in a combustion chamber (Case 3).

Figure 20 :

 20 Figure 20: Relative standard deviation as a function of the number of realizations N for three characteristic points of the channel flow configuration (Case 1): A (red), B (black) and C (blue). Circles: MC simulations. Diamonds: RQMC simulations. Grey lines: convergence rate of MC (dashed line) and RQMC (dotted line) results.

Figure 21 :

 21 Figure 21: Relative standard deviation as a function of the number of realizations N for two characteristic points of the sooting jet flame (Case 2): D (red) and E(black). Circles: MC simulations. Diamonds: RQMC simulations. Grey lines: convergence rate of MC (dashed line) and RQMC (dotted line) results.

Figure 20 Figure 22 :

 2022 Figure20shows that the MC standard deviation (circles) is larger than the RQMC one (diamonds) for the high and low temperature probes (A and C) for every number of rays N considered. Instead,

Figure 23 :

 23 Figure 23: Field of the required number of rays in controlled-convergence Test 2 with MC (top) and RQMC (bottom) for the channel flow configuration (Case 1).

Figure 24 :

 24 Figure 24: Field of the required number of rays in controlled-convergence Test 2 with MC (top) and RQMC (bottom) in a sooted jet flame application (Case 2).

Figure 25 :

 25 Figure 25: Field of the required number of rays in controlled-convergence Test 2 with MC (top) and RQMC (bottom) for a combustion chamber (Case 3).

590 σ 2 i

 2 [Q(N )] computed from Eqs. (11) for MC and (14) for RQMC, and the local computational time, T CP U i

Figure 26 :

 26 Figure 26: 2D map of the ratio between the efficiency of RQMC and the efficiency of MC methods for a channel flow application (top), a jet flame (middle) and a confined combustion chamber (bottom).

3 . 2 ,

 32 consists in sampling the frequency according to the emission at the highest temperature of the computational domain. Both ERM and OERM are combined to MC and RQMC formalisms and the four corresponding combinations, MC-ERM, MC-OERM, RQMC-ERM and RQMC-OERM, are compared in the 1D test case studied in Sec. 3.4, called Case b.The improvement of the convergence in cold regions is outlined by comparing the methods at a point close to the wall at the temperature T = 550 K. Figure27displays the relative standard deviation as a function of the number of realizations achieved by MC (circles) and RQMC (diamonds) when combined to ERM (blue) and OERM (red). Grey dashed lines highlight the convergence rate of MC: it can be seen that when MC is combined to OERM the line of slope -1/2 is shifted down because of the lower standard deviation obtained with OERM. Grey and black dotted lines instead represent the convergence rate of RQMC when combined with ERM and OERM, respectively. It is worth noticing the improvement accomplished by RQMC and OERM: the combination of RQMC and importance sampling technique shifts down the relative standard deviation obtained from RQMC-ERM, as also seen in MC results (dashed grey line), leading to the lowest standard deviation achievable among the four investigated combinations. The convergence rate of RQMC can then be indeed further accelerated when RQMC is combined to the OERM method.

Figure 27 :

 27 Figure 27: Relative standard deviation as a function of the number of realizations N for a cold point close to the wall in the slab configuration achieved in MC (circles) and RQMC (diamonds) combined with ERM (blue) and OERM (red). Grey dashed line: convergence rate of MC. Grey dotted line: convergence rate of RQMC-ERM. Black dotted line:

  and HITEMP 2010 for the ones of H 2 O [49]. The whole considered spectrum ranges from 150 to 9200 cm -1 . 44 spectral bands are considered for H 2 O, whose width varies from 50 to 400 cm -1 . CO 2 spectrum overlaps the H 2 O one and CO 2 absorbs in only 17 of the 44 bands. This leads to a total of 1022 quadrature points: 7 2 quadrature points for each of the 17 bands where both species absorb plus 7 quadrature points for the remaining 27 bands where only H 2 O participates.

Table 1 :

 1 Numerical set-up of the tests with a fixed number of realizations, Test 1, and at controlled convergence, Test 2, for the three applications.

		Test 1	Test 2
	Case Number of Rays Rel. Std	Abs. Std
	1	5 000	1%	3% of P R max
	2	5 000	1%	1% of P R max
	3	5 000	3%	3% of P R max

.4. Numerical set-up for MC and QMC simulations For each application, radiative heat transfer simulations are performed using both MC and RQMC methods. For each method two different tests are carried out. The first test, Test 1, consists in

Table 2 :

 2 Wall clock time in seconds needed for MC (3 rd column) and RQMC (4 th column) simulations (Test 2) for the retained configurations using 168 cores. The number of mesh nodes in each case is indicated in the second column

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 643134. This work was also supported by the Air Liquide, CentraleSupélec and CNRS Chair on oxy-combustion and heat transfer for energy and environment and by the OXYTEC project, grant ANR-12-CHIN-0001 of the French Agence Nationale de la Recherche. The numerical study was carried out thanks to the HPC resources of CINES under the allocations 2016-020164 and 2018-A0042B10159 made by GENCI, and the HPC resources from the "Mésocentre" computing center of CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and grants from Ministry of Research via the Contrat Plan État-Région.