Optimized Emission-based Reciprocity Monte Carlo Method to speed up computation in complex systems

Keywords: ]. β Monte Carlo, frequency distribution function, Reciprocity Method, radiative power, statistic error

If radiation plays an important role in many engineering applications, especially in those including combustion systems, influence of radiation on turbulent flows, particularly on the turbulent boundary layers, is still not well known.The objective is here to perform a detailed study of radiation effect on turbulent flows. An optimized emission-based reciprocal (OERM) approach of the Monte-Carlo method is proposed for radiation simulation using the CK model for radiative gas properties. OERM allows the uncertainty of results to be locally controlled while it overcomes the drawback of the original emission-based reciprocity approach by introducing a new frequency distribution function that is based on the maximum temperature of the domain. Direct Numerical Simulation (DNS) has been performed for turbulent channel flows under different pressure, wall temperatures and wall emissivity conditions. Flow field DNS simulations are fully coupled with radiation simulation using the OERM approach. The role of radiation on the mean temperature field and fluctuation field are analyzed in details. Modification of the mean temperature profile leads to changes in wall conductive heat fluxes and new wall laws for temperature when radiation is accounted for. The influence on temperature fluctuations and the turbulent heat flux is investigated through their respective transport equations whose balance is modified by radiation. A new wall-scaling based on the energy balance is proposed to improve collapsing of wall-normal turbulent flux profiles among different channel flows with/without considering radiation transfer. This scaling enables a new turbulent Prandtl number model to be introduced to take into account the effects of radiation. In order to consider the influence of radiation in the near-wall region and predict the modified wall law, a one-dimensional wall model for Large Eddy Simulation (LES) is proposed. The 1D turbulent equilibrium boundary layer equations are solved on an embedded grid in the inner layer. The obtained wall friction stress and wall conductive flux are then fed back to the LES solver. The radiative power term in the energy equation of the 1D wall model is computed from an analytical model. The proposed wall model is validated by a comparison with the former DNS/Monte-Carlo results. Finally, two criteria are proposed and validated. The first one is aimed to predict the importance of wall radiative heat flux while the other one predicts whether a wall model accounting for radiation in the near wall region is neces-

I would like to thank my supervisors, Prof. Olivier GICQUEL and Prof. Jean TAINE, for their support during my Ph.D. study. I could not complete my thesis without their enlightening instruction, impressive kindness and patience. I am also grateful to my co-supervisor, Ronan VICQUELIN, for his patient advice, warm help and edits of the documents. I was very lucky to be his first Ph.D. student. I would also thank all my friends in the laboratory for their help in my study and happy time together. Among them, Dr. Jingxuan LI, who was also my college when I was in China, has gave me a lot of help in both my study and life. He was the kind person who went to the airport and picked me up when I arrived in France three years ago. Moreover, I want to thank Dr. Fabien THOLIN for his warm help. I will miss the moment when I played badminton with him. And I am also appreciated for the friendship of Jean-francois BOUR-GOUIN who was in the same office as me. He is a very funny and kind person and he has brought me a lot of laugh in the past three years. I want to acknowledge Dr. Philippe Riviere for providing the database of the Correlated-K model for my simulation and I am grateful to Dr. Vincent Moureau and the SUCCESS scientific group for providing the code YALES2 and support. I also express thanks to all the other professors, researchers and secretaries who have helped in the past three years.

I am thankful to all my Chinese friends in ECP (Ecole Centrale Paris) who encouraged and supported me in these three years. Because of your company, I didn't feel lonely even during some traditional Chinese festivals when we are supposed to be with families.

I also want to thank the China Scholarship Council (CSC) for granting me a doctoral fellowship so that I can have this wonderful experience in France.

I am sincerely grateful to my parents, Zhouhua ZHANG and Liqun HUANG, my little sister, Yuling ZHANG and all the other family members. Thank you for consistently encouraging me to be brave, independent and optimistic, and Finally but not least, I want to thank boyfriend, Ying ZHU. Although there was a long physical distance between us (he was in China while I was in France), he was always there cheering me up and listening to me during all the difficult moments.

Résumé

Le rayonnement joue un rôle important dans de nombreuses applications industrielles, en particulier celles mettant en jeu un processus de combustion. Cependant, son influence sur les écoulements turbulents, notamment les couches limites, n'est pas encore bien connu. L'objectif est ici d'analyser en détail l'effet du rayonnement sur les écoulements turbulents. Utilisant le modèle CK pour décrire les propriétés radiatives des gaz, une approche optimisée de la méthode de Monte-Carlo (OERM) basée sur l'émission et le principe de réciprocité est développée. La méthode OERM permet de contrôler localement l'imprécision des résultats tout en corrigeant l'inconvénient de la méthode originale en introduisant une nouvelle fonction de répartition des fréquences basée sur la température maximale du domaine. Plusieurs écoulements de canal plan turbulent sous différentes conditions de pression, de températures et d'émissivités de parois sont résolus par simulation numérique directe (DNS). Les simulations DNS de l'écoulement et du champ de rayonnement par la méthode OERM sont entièrement couplées. L'impact du rayonnement sur le champ moyen de température et ses fluctuations est analysé en détail. La modification du profil de température moyenne induit un changement des flux de chaleur conductifs aux parois et de nouvelles lois de paroi pour la température lorsque le rayonnement est pris en compte. L'impact sur les fluctuations de température et le flux de chaleur par transport turbulent est étudié au travers de leurs équations de transport respectives dont l'équilibre est modifié par le rayonnement. Une nouvelle normalisation (wall-scaling) basée sur le bilan d'énergie est proposée pour améliorer le recouvrement des profils obtenus sous les différentes configurations étudiées avec et sans transfert radiatif. Cette normalisation permet d'introduire un modèle pour le nombre de Prandtl turbulent lorsque le rayonnement est pris en compte. Afin de prédire l'effet du rayonnement sur la zone proche paroi et sa modification des lois de paroi, un modèle de paroi pour la simulation aux grandes échelles (LES) est développé. Les équations 1D de couche limite turbulente à l'équilibre sont résolues sur une grille intégrée au maillage LES pour traiter la zone interne. La contrainte pariétale et le flux de chaleur conductif obtenus sont renvoyés au code LES. La puissance radiative dans l'équation d'énergie du modèle de paroi 1D est évaluée à partir d'un modèle analytique. Le modèle de paroi est validé par comparaison avec les résultats des calculs couplés Deux critères sont finalement proposés et validés. Le premier a pour but de prédire l'importance du flux radiatif pariétal tandis que le second détermine si un modèle de paroi prenant en compte l'effet du rayonnement dans la zone interne de la couche limite est nécessaire. Une étude paramétrique est ensuite réalisée où un modèle k-ǫ et un modèle de nombre de Prandtl turbulent sont utilisés pour estimer les champs moyens de vitesse et température d'écoulements de canal plan sous différentes conditions. Les valeurs des critères obtenues sont analysées puis comparées.

Nomenclature

Roman Symbols

Introduction

Thermal radiation is one of the most important transfer modes. Any object emits thermal radiation. As temperature increases, radiative transfer rate increases markedly. Thereby, it plays an extremely important role in many engineering applications, especially in combustion systems and in plasmas with high operating temperature and pressure conditions.

Influence of radiation transfer

In most of industrial furnaces, due to a large combustor size, radiation often plays a predominant role over convective heat transfer. Fig. 1 shows a detailed energy flow diagram for an industrial furnace. Among the 35% useful output, only 4% is contributed by convection while all the other heat transfer is from radiation. Therefore, prediction and control of radiative transfer is crucial for furnace performance.

In modern gas turbine combustors, as shown in Fig. 2, the peak temperatures of the hot gases in primary zone can be higher than 2100K, whereas for the materials commonly used for liner and dome, the maximum operating temperature should not exceed 1100K [START_REF] Lefebvre | Gas Turbine Combustion[END_REF]. Consequently, protecting liner walls and domes from hot gases is critical and challenging for combustor designers.

Generally, film-cooling air is applied to protect combustor walls. To control the quantity of cooling air, information about the heat fluxes on walls is necessary.

In the presence of film cooling, a sizable fraction of the heat transferred from gas to combustor walls is by radiation. And this importance of radiative flux is even higher in modern gas turbine as operating pressure and temperature increase. Moreover, in order to release more air for mixing and controlling of pollutant emission, the amount of air employed for cooling of liner walls must be minimized. The determination of the minimum quantity requires a sound knowledge of radiation.

The pollutant emission of combustors is also influenced by radiation due to its sensitivity to local temperature. It follows that accurate assessment of pollutant emission necessitates precise prediction of local radiative transfer rate.

Coupled turbulent convective transfer and radiative transfer

In most of the practical systems involving heat transfer process, the flow is turbulent and hence radiation is coupled with turbulent convection. For in- stance, as shown in Fig. 3, the basic heat transfer processes for a liner wall of combustor include: heating by radiation and turbulent convection from the hot combustion gas, cooling by radiation to the outer casing and by turbulent convection to the annulus air and heat transferring through the wall by conduction. The liner-wall temperature is then determined by the balance of these processes. Among these heat transfer processes, radiative transfer differs from conduction and convection. The latter two are short-range phenomena and commonly the energy conservation could be performed in an "infinitesimal volume". By contrast, radiation is a long-range phenomenon and an integro-differential equation is required in order to represent it. The difficulty of radiation problem is even further enhanced by the need to consider radiative properties of participative gases, which generally vary significantly with wave number. Generally, the radiation problem can be investigated experimentally or numerically. However, due to formidable difficulties encountered, radiation has not been subjected to extensive experimental investigations, especially in combustion at high pressure. Fortunately, numerical simulation appears as an extremely useful tool.

Numerical simulation of radiation and turbulence

In early stage of combustor design, global radiation models were widely applied to assess radiation contribution. The simple, global model based on mean beam length (Hottel hemisphere) concept does not consider the non-uniformities of temperature and species concentration, and hence could lead to inaccurate results.

Thanks to an increasing understanding of radiative transfer, various methods have then been developed, including zonal method, ray tracing method, direct ordinate model, spherical harmonics model, Monte-Carlo method, etc.. Among them, Monte-Carlo is one of the most promising methods since it can be easily employed in any complicated system and can incorporate all the important radiative effects without any assumption and large increment of computational effort. However, its computational cost is high since it needs a large number of rays to obtain statistically meaningful results.

Regarding the simulation of turbulence, a proper treatment of near wall region is alway an important issue, especially in Large Eddy Simulation (LES). As shown in Fig. 4 (a), for resolving the near wall region, a fine grid is required in the wall vicinity and this yields a high resolution requirement. To alleviate this constraint, wall models are employed to treat the near wall region and a coarse grid can then be applied as shown in Fig. 4 (b). The chosen wall model uses the velocity u and the temperature T at LES grid point y as outer boundary conditions and the wall friction stress τ w and the wall conductive heat flux q w are then calculated and sent back to LES solver. Wall-modeled LES is more and more applied in engineering applications characterized by a high Reynolds number.

Various numerical investigations have also been performed in order to investigate the coupling between turbulent flows and radiation. It has been found that, in reactive flows, radiation can affect the mean or local temperature profile [START_REF] Dos Santos | Coupled large eddy simulations of turbulent combustion and radiative heat transfer[END_REF][START_REF] Damien | Analysis of the interaction between turbulent combustion and thermal radiation using unsteady coupled LES/DOM simulations[END_REF][START_REF] Liu | Effects of gas and soot radiation on soot formation in counterflow ethylene diffusion flames[END_REF], NO x emission [START_REF] Barlow | Scalar profiles and no formation in laminar opposed-flow partially premixed methane/air flames[END_REF]) and soot volume fraction [START_REF] Liu | Effects of gas and soot radiation on soot formation in counterflow ethylene diffusion flames[END_REF]. For example, in the study of dos [START_REF] Dos Santos | Coupled large eddy simulations of turbulent combustion and radiative heat transfer[END_REF], LES of a turbulent premixed flame has been performed and a ray-tracing method has been used for the calculation of the radiation field. They have reported that radiation decreases the maximum temperature by about 200 K (as shown in Fig. 5), homogenizes the mean temperature field and strongly enhances the flame dynamics. [START_REF] Damien | Analysis of the interaction between turbulent combustion and thermal radiation using unsteady coupled LES/DOM simulations[END_REF] have also found that radiation decreases the mean temperature and changes the local extreme temperature in a laboratory flame configuration, while [START_REF] Liu | Effects of gas and soot radiation on soot formation in counterflow ethylene diffusion flames[END_REF] have observed a change of flame peak temperature and soot volume fraction due to radiation in counterflow ethylene diffusion flames. Meanwhile, regarding the effect of radiation on non-reactive flows, [START_REF] Soufiani | Radiation effects on turbulent heat transfer in channel flows of infrared active gases[END_REF] have found that, in turbulent channel flow, the temperature profile and the conductive heat flux on walls are strongly affected by radiation. [START_REF] Gupta | Large-Eddy Simulation of Turbulence-Radiation Interactions in a Turbulent Planar Channel Flow[END_REF] have also observed a similar effect of radiation on the mean temperature profile. However, in the study of [START_REF] Amaya | Unsteady coupling of Navier-Stokes and radiative heat transfer solvers applied to an anisothermal multicomponent turbulent channel flow[END_REF], a weak effect of radiation is found in a similar channel flow, which can be attributed to the small optical thickness of the medium. Furthermore, great attention has been given to the interaction between turbulence and radiation (TRI). Experimental studies [START_REF] Gore | Spectral and total radiation properties of turbulent carbon monoxide/air diffusion flames[END_REF][START_REF] Ji | Thermal radiation properties of turbulent lean premixed methane air flames[END_REF][START_REF] Zheng | Measurements and inverse calculations of spectral radiation intensities of a turbulent ethylene/air jet flame[END_REF], theoretical analysis [START_REF] Burns | Turbulence radiation interaction modeling in hydrocarbon pool fire simulations[END_REF][START_REF] Baum | Radiation and velocity fields induced by localized temperature fluctuations[END_REF][START_REF] Soufiani | Radiation effects on turbulent heat transfer in channel flows of infrared active gases[END_REF]) and numerical simulations [START_REF] Wu | Direct numerical simulation of turbulence/radiation interaction in premixed combustion systems[END_REF][START_REF] Deshmukh | Direct numerical simulation of turbulence-radiation interactions in homogeneous nonpremixed combustion systems[END_REF][START_REF] Deshmukh | Direct numerical simulation of turbulence-radiation interactions in a statistically onedimensional nonpremixed system[END_REF][START_REF] Dos Santos | Coupled large eddy simulations of turbulent combustion and radiative heat transfer[END_REF]) have been carried out to investigate the effect of TRI in different systems and comprehensive reviews are available in [START_REF] Coelho | Numerical simulation of the interaction between turbulence and radiation in reactive flows[END_REF][START_REF] Coelho | Turbulence-Radiation Interaction: From Theory to Application in Numerical Simulations[END_REF]. TRI arises from the high non-linear coupling between fluctuation of temperature and species concentration and fluctuation of radiation. It consists of two parts: effect of turbulence on radiation field and vice-versa.

Regarding the effect of turbulence on radiation, it is observed that turbulence leads to an increase in the medium transmissivity [START_REF] Jeng | Radiative heat fluxes near turbulent buoyant methane diffusion flames[END_REF][START_REF] Gore | Spectral and total radiation properties of turbulent carbon monoxide/air diffusion flames[END_REF], the radiative power [START_REF] Coelho | Detailed numerical simulation of radiative transfer in a nonluminous turbulent jet diffusion flame[END_REF][START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF]) and radiative energy loss [START_REF] Li | Importance of turbulence radiation interactions in turbulent diffusion jet flames[END_REF][START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF]). Coelho et al. have reported that, in a non-luminous turbulent jet diffusion flame (Coelho et al. 2003), TRI enhances the radiative energy loss by a factor of 30% while a similar change in a luminous turbulent flame is also revealed in Ref. [START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF]).

Introduction

Moreover, the contribution of TRI on the mean radiation field is constituted with different contributions: the temperature self-correlation, the absorption coefficient self-correlation, the absorption coefficient-temperature correlation and the absorption coefficient-radiative intensity correlation [START_REF] Coelho | Numerical simulation of the interaction between turbulence and radiation in reactive flows[END_REF].

In order to isolate and quantify individual contributions to TRI in a statistically one-dimensional premixed combustion system, [START_REF] Wu | Direct numerical simulation of turbulence/radiation interaction in premixed combustion systems[END_REF] have performed a direct numerical simulation (DNS) coupled with a Monte-Carlo method. Their results have revealed that the temperature self-correlation contribution is only dominant in the case with smallest optical thickness in their study. Among the three correlation terms of the absorption coefficient, the ones with the Planck function and with the incident radiative intensity are not negligible, even in the most optically thin case, while for a case with intermediate value of optical thickness, the three correlations are all significant. Deshmukh et al. [START_REF] Deshmukh | Direct numerical simulation of turbulence-radiation interactions in homogeneous nonpremixed combustion systems[END_REF][START_REF] Deshmukh | Direct numerical simulation of turbulence-radiation interactions in a statistically onedimensional nonpremixed system[END_REF] have also studied different contributions to TRI with DNS in a statistically homogeneous isotropic nonpremixed combustion system and a one-dimensional turbulent non-premixed flame. Only the latter one was fully coupled with radiation. By contrast to the former described effects, only a few studies have been devoted to the effect of radiation on turbulence. Among them, [START_REF] Soufiani | Temperature turbulence spectrum for high-temperature radiating gases[END_REF] have carried out a theoretical analysis of the influence of radiation on thermal turbulence spectra and it has been concluded that radiation acts as a dissipation term and it could smooth the intensity of temperature fluctuation and modify the structure of temperature spectra. [START_REF] Damien | Analysis of the interaction between turbulent combustion and thermal radiation using unsteady coupled LES/DOM simulations[END_REF] have also reported that radiation modifies the temperature fluctuation while the energy and frequency distribution is homogenized. Moreover, [START_REF] Ghosh | Effects of radiative heat transfer on the structure of turbulent supersonic channel flow[END_REF] have presented that radiation modifies the Reynolds stress and affects the turbulence structure in supersonic shear layers in a non-negligible way.

The objective of this thesis

Despite of the various studies of radiation effect on turbulent flows, there still exist several fields of research:

• No general trend or physical understanding of the influence of radiation in turbulent boundary layer is obtained;

• The existing wall models used in LES do not consider radiation effect which could be very important in some conditions;

• There is no adequate principle to assess the necessity of a radiation simulation under certain conditions in practical applications.

Introduction

23

The main objective of this thesis is to explore better understanding and prediction of the influence of radiation in turbulent flows. The structure of the manuscript is detailed in the following.

Structure of this manuscript

Chapter 1 Notions related to radiation and various approaches to represent participating radiative properties and to solve the radiative transfer equation are introduced. Different Monte-Carlo methods, especially the Reciprocal Monte-Carlo method, are detailed;

Chapter 2 An Optimized Emission-based Reciprocity Monte Carlo Method (OERM) has been developed and validated by comparison, in benchmark cases, with analytical reference results. In this method, the frequency distribution function, generally equal to the emission distribution function at the emitting cell temperature, has been replaced by a distribution function associated with the maximum temperature within the domain. The real emission distribution function of a cell at any temperature is then obtained by applying a corrective factor to each shot;

Chapter 3 After an introduction of the balance equations of turbulent flows and of the three numerical approaches (DNS, LES and RANS), wall functions associated with the boundary layers are presented. The code YALES2 used for flow field simulation and the coupling between this code and the Monte-Carlo method are also detailed;

Chapter 4 Direct numerical simulations of a turbulent flow in a two-dimensional channel are performed. The fluid solver is coupled to the reciprocal Monte Carlo method OERM when radiation is taken into account. Effects of gas-gas radiation and gas-wall radiation on the mean temperature and flux field are studied separately. The dependence of radiation effects on the set of conditions (temperature level, wall emissivity, pressure, Reynolds number) is also analyzed;

Chapter 5 This chapter focuses on radiation effects on fluctuation field in turbulent channel flow. The influence of radiation on enthalpy root-mean-square, on turbulent heat flux and on the budget of the corresponding transport equation are analyzed. To improve the agreement of results of non-dimensional turbulent quantities between flows without and with radiation, a new radiationbased scaling is proposed. The influence of radiation on turbulent Prandtl number is also presented and a model based on the new scaling is proposed;

Chapter 6 In order to take into account the effect of radiation on the temperature wall law in turbulent boundary layers, a new wall model for large eddy Introduction simulation (LES) is proposed. After a separate validation of different components of the proposed wall model, fully coupled wall-modeled LES/Monte-Carlo simulations of turbulent channel flow are performed. The obtained results are compared with DNS/Monte-Carlo results from chapter 4 to assess the accuracy of this new model;

Chapter 7 A criterion regarding the importance of radiative heat flux on walls is proposed to assess the necessity of performing radiation simulation under a given flow condition. Similarly, another criterion based on the change of non-dimensional temperature (scaled in wall units) due to radiation, is also developed to predict whether the wall model accounting for radiation is required for a given flow. The validity of these two criteria is confirmed from an analysis of fully-coupled DNS/Monte-Carlo results of several turbulent channel flows. Mean flow fields of channel flows under different flow conditions are then resolved with a k-ǫ model and a model for the turbulent Prandtl number. The obtained criterion results are presented and analyzed.

Introduction to radiation simulation

In this section, the radiative transfer equation (RTE) is detailed and different approaches to represent participative radiative properties and to solve the RTE are discussed. Finally, a detailed description of different Monte-Carlo methods, especially the Reciprocal Monte-Carlo method, is also presented.

Radiative transfer

Thermal radiative transfer is commonly described as the heat transfer caused by electromagnetic waves. As one of the three heat transfer modes, radiation transfer is quite different from the other two modes: conduction and convection. As we know, the latter two can be both described by using a partial differential equation since they are short-range phenomena and commonly the energy conservation could be performed in an "infinitesimal volume". By contrast, radiation is a long-range phenomenon and the energy conservation need to be carried out over the entire enclosure under consideration. To represent the radiation, an integro-differential equation is retained with seven independent variables: three spatial coordinates, two angle coordinates, one spectral variable and the time. Although generally the time dependence is omitted in most practical applications since the transfer speed of radiation beam is superior than the typical velocity of such flows by several order of magnitude, there remains six independent variables. The radiation problem is further complicated by the need to calculate radiative properties of the medium. The absorption coefficient of radiant gases, especially CO 2 and H 2 O, is constituted by hundreds of thousands of spectral lines and it varies significantly with wave number. This makes the radiation a difficult problem to resolve. On the other hand, radiative transfer rates are crudely proportional to differences of temperature to the fourth (or higher) power. Therefore, as temperature increases, radiative transfer becomes more important and may be totally dominant over conduction and convection at very high temperatures. As a consequence, in combustion applications where the temperature is high, radiation is important and needs to be modeled accurately.

Radiative Transfer equation

Radiative intensity is one of the most used quantities to describe the radiative phenomenon. It is defined as the radiative energy flux per unit solid angle and per surface area normal to the rays. The total radiative intensity is the intensity over the entire spectrum while the spectral radiative intensity is the one at given frequency and per unit frequency interval.

As shown in Fig. 1.1, in an elemental cylindrical control volume which contains emitting-absorbing-scattering medium, the incident spectral radiative intensity I ν , with a direction s and a wave-number ν, could be attenuated by absorption of the medium in the control volume and by the deviation of the beam from the original direction s (out-scattering). On the other hand, the emission of the medium and in-scattering from the radiative intensity in other propagation directions s i to the direction s contribute to amplify the radiative intensity.

From an analysis of energy balance of the control volume in the direction s the radiative transfer equation (RTE), for homogeneous and isotropic participative media with a refractive index equal to 1 can be written as

dI ν ds = -κ ν I ν -σ sν I ν + κ ν I 0 ν + σ sν 4π 4π I ν (s i )Φ ν (s, s i )dΩ i , (1.1)
where κ ν and σ sν are the spectral absorption and scattering coefficients respectively and Φ i the so-called scattering phase function which represents the probability that a beam from direction s i , that is scattered within is scattered to direction s within dV (as shown in Fig. 1.2). Ω denote the solid angle. The four terms on the right side of Eq. (1.1) denote the extinction by absorption, by out-scattering, augmentation by emission and by in-scattering respectively. The dependence of intensity on time is neglected in this equation since commonly the speed of photons exceeds velocities of other terms by many order of magnitude. Introducing the extinction coefficient β ν , the single scattering albedo ω ν and the non-dimensional optical thickness τ ν as

β ν = κ ν + σ sν , (1.2
) where the source function S ν (τ ν , s) is defined as

ω ν = σ sν κ ν + σ sν , (1.3)
S ν (τ ν , s) = (1 -ω ν )I 0 ν + ω ν 4π 4π I ν (s i )Φ ν (s, s i )dΩ i .
(1.6)

An integration of Eq. (1.5) from s ′ = 0 to s ′ = s gives rise to the integral formulation of the RTE as

I ν (τ ν ) = I ν (0)e -τν + τν 0 S ν (τ ′ ν , s)e -(τν -τ ′ ν ) dτ ′ ν .
(1.7)

Here the first term on the right side represents the incident radiative intensity at s ′ = 0 attenuated along the trajectory of the beam until s ′ = s while the second term is the contribution of the source term (by emission or scattering) of every element and attenuated over the distance τ ν -τ ′ ν . If the medium is non-scattering (σ sν = 0), the differential and integral formulations reduce to the following forms respectively: (1.9)

dI ν ds = κ ν (I 0 ν -I ν ), (

Boundary condition

Generally, the radiative intensity leaving a wall surface which surrounds the participation medium, can be specified and employed as boundary conditions for the RTE.

For diffusely emitting and reflecting opaque surfaces, the exiting radiative intensity at position r w on the surface is independent of the direction, and it can be expressed as the sum of emitted intensity and reflected intensity:

I ν (r w , s) = I ν (r w ) = ǫ(r w )I 0 ν (r w ) + ρ(r w )

H ν (r w ) π , (1.10)
where H ν is the hemispherical irradiation defined as (1.11) where, as indicated in Fig. 1.3, n w is the surface normal unit vector and s ′ is the direction unit vector of the incoming intensity I ν (r w , s ′ ).

H ν (r w ) = nw•s ′ <0 I ν (r w , s ′ )|n w • s ′ |dΩ ′ ,
It is often considered that the radiative properties are a combination of diffuse and specular ones, and then Eq. (1.10) changes to:

I ν (r w , s) = ǫ(r w , s s )I 0 ν (r w ) + ρ d (r w )
H ν (r w ) π + ρ s (r w , s s )I ν (r w , s s ), (1.12)

where ρ s and ρ d are the specular and diffuse components of the reflectance respectively. The wall emissivity is ǫ = 1 -ρ, where the total reflectance ρ, while not diffuse, can be adequately represented by a combination of a diffuse and a specular component as

ρ = ρ s + ρ d .
(1.13)

Besides, the direction s s and s have a relationship as

s s = s -2(s • n w )n w . ( 1

Radiative source term

In order to correlate radiation field with flow field, a source term P R (W/m 3 ), named radiative power per unit volume, is added in the energy balance equation of the flow. The radiative power per unit volume P R is defined as the negative derivative of radiative flux q R :

P R = -▽ •q R = ▽ • +∞ 0 q R ν dν (1.15)
where q R ν , the spectral radiative flux crossing a surface with normal unit vector n, is expressed as

q R ν = q R ν • n = 4π I ν n • sdΩ (1.16)
And the divergence of the spectral radiative flux q R ν writes

▽ • q R ν = ▽ • 4π I ν s dΩ (1.17)
Integrating Eq. (1.1) over the entire solid angle gives rise to

▽ • q R ν = ▽ • 4π I ν s dΩ = 4πκ ν I 0 ν - 4π β ν I ν dΩ + σ sν 4π 4π I ν (s i ) 4π Φ ν (s, s i )dΩ dΩ i (1.18)
In case of a non-scattering isotropic medium, this equation can be simplified as:

▽ • q R ν = κ ν (4πI 0 ν - 4π I ν dΩ) (1.19)
This expression implies that the variation of radiative energy inside a control volume is the difference between emitted energy and absorbed energy.

Radiative properties of participative gas

The prediction of radiative properties of medium is one of the most important factors that make describing the radiation phenomenon a difficult task. Participative gases, particularly CO 2 and H 2 O which are typical products of combustion system, emit and absorb radiative energy selectively. Consequently, it is important to take into account the wave-number dependence of radiative properties. Exact prediction results can only be obtained with line-by-line calculations [START_REF] Taine | A line-by-line calculation of low-resolution radiative properties of CO 2 -CO-transparent nonisothermal gases mixtures up to 3000 k[END_REF][START_REF] Hartmann | Line-by-line and narrow-band statistical model calculations for h2o[END_REF] which represents each discrete absorption-emission line of the infinite spectrum by its strength and half-width. However, due to its requirement of vast amounts of computer resources, Line-byline method is only used as benchmarks for the validation of more approximate 1.1. RADIATIVE TRANSFER spectral models [START_REF] Modest | radiative heat transfer[END_REF]). In practical applications, it is not feasible to account for each individual spectral line and hence a model for the spectral calculation is needed. In the literatures, there exist three different types of approaches for the representation of radiative properties: Narrow-band models, wide-band models and Global approach.

• Narrow-band models Over a narrow spectrum band, although the absorption coefficient varies rapidly with wave number, other radiative quantities such as the blackbody intensity, doesn't change appropriately. Hence, it is reasonable to replace the actual absorption coefficient by smoothened values appropriately averaged over a narrow spectral range. In order to calculate the averaged value, some information on the spacing of individual lines within the narrow band and on their relative strengths is needed. Various narrow band models have been proposed for this purpose. In the Elsasser model [START_REF] Elsasser | Heat transfer by infrared radiation in the atmosphere[END_REF], it is assumed that the lines are of uniform intensity and are equally spaced while the Goody or statistical model [START_REF] Goody | A statistical model for water-vapour absorption[END_REF] postulates a random exponential line intensity distribution and a random line position selected from a uniform probability distribution. Both of these models yield reasonably accurate results for practical engineering calculations.

An alternative to the "traditional" narrow band models is the so-called "correlated κ-distribution" (CK model) [START_REF] Goody | The correlated-k method for radiation calculations in nonhomogeneous atmospheres[END_REF]. Provided that in the narrow spectral range, the averaged band quantities are not affected by the precise location of the frequency, the absorption coefficient could be reordered, resulting in a smooth monotonically increase function of absorption coefficient vs. artificial wavenumber. This function provides an efficient and accurate scheme for the spectral integration of radiative quantities of interest. A more detailed description of this model could found in Ref. [START_REF] Soufiani | High temperature gas radiative property parameters of statistical narrow-band model for H 2 O, CO 2 and CO, and correlated-K model for H 2 O and CO 2[END_REF][START_REF] Taine | Gas ir radiative properties: From spectroscopic data to approximate models[END_REF]). The CK model has been widely used in radiation studies [START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF][START_REF] Zhang | Monte Carlo References method of radiative transfer applied to a turbulent flame modeling with LES[END_REF]) and is also used in this thesis for the representation of the gas radiative properties.

• Wide-band models In the 70s, a correlated approach at scales of a spectroscopic wide band has been developed by [START_REF] Edwards | Comparison of models for correlation of total band absorption[END_REF][START_REF] Edwards | Molecular gas band radiation[END_REF]). Due to the fact that the necessary calculations are relatively simple, this wide band model was very popular in the past. However, nowadays, it is not commonly used anymore because of its low correlational accuracy.

• Global models Among various global models, the models based on the distribution of the [START_REF] Taine | Gas ir radiative properties: From spectroscopic data to approximate models[END_REF]. These models are commonly called the weighted-sum-of-gray-gases (WSGG) models. The original version of the WSGG model is proposed by [START_REF] Hottel | Radiative transfer[END_REF], where the total emissivity of the medium is expressed as the summation of several virtual gray gases with weighting factors depending on the temperature. The main interest of WSGG model is its low computational cost. However, their application is generally limited to problems with grey walls and/or particles.

More accurate global model could be found in Refs. [START_REF] Denison | A Spectral Line-Based Weighted-Sum-of-Gray-Gases Model for Arbitrary RTE Solvers[END_REF] and [START_REF] Pierrot | A fictitious-gas-based absorption distribution function global model for radiative transfer in hot gases[END_REF] where the spectral line-based weightedsum-of-grey-gases model (SLW) and the absorption distribution function model with fictitious gases (ADFFG) are proposed respectively.

Numerical approaches to resolve RTE

In order to solve the RTE, various approximated models have been proposed in the literature. Depending on the form of the resolved RTE, these models can be divided into two classes: differential and integral methods.

In the catalogue of differential methods, the differential form of RTE is resolved by using an approximation of the incident intensity. It includes the discrete ordinate model [START_REF] Fiveland | Discrete-ordinates solutions of the radiative transport equation for rectangular enclosures[END_REF], finite volume model [START_REF] Raithby | A Finite-Volume Method for Predicting a Radiant Heat Transfer in Enclosures With Participating Media[END_REF], spherical harmonics model [START_REF] Mengüç | Radiative Transfer in Axisymmetric, Finite Cylindrical Enclosures[END_REF]) and discrete transfer methods [START_REF] Lockwood | A new radiation solution method for incorporation in general combustion prediction procedures[END_REF].

On the other hand, integral methods resolve the integral formulation of RTE and the mostly used methods include the zonal method [START_REF] Hottel | Radiative transfer[END_REF], ray tracing method [START_REF] Iacona | Radiation transfer modeling by a ray-tracing method and a ck or k approach[END_REF]Coelho et al. 2003), boundary element [START_REF] Bialecki | Solving heat radiation problems using the boundary element method[END_REF] and Monte-Carlo method [START_REF] Howell | The Monte Carlo method in radiative heat transfer[END_REF]. A detailed description about these methods can be found in some comprehensive references or books about radiative transfer, such as Ref. [START_REF] Modest | radiative heat transfer[END_REF][START_REF] Viskanta | Radiation heat transfer in combustion systems[END_REF]. Among these models, Monte-Carlo method is a statistical method, where the history of a large number of photon bundles is traced as they travel through the enclosure and the radiative quantities are statistically calculated. This model is well-known for its high accuracy of the solution. In this thesis, the Monte-Carlo method is retained in order to predict the result with a good accuracy.

A detailed description of Monte-Carlo methods will be presented in the next paragraph.

Monte-Carlo Method

The Monte-Carlo method (MCM), based on the statistical characteristics of physical processes, has been extensively used in engineering and science investigations. The application of MCM in thermal radiative transfer was pioneered by [START_REF] Fleck | The calculation of nonlinear radiation transport by a monte carlo method: Statistical physics[END_REF], [START_REF] Howell | Monte carlo solution of thermal transfer through radiant media between gray walls[END_REF] who studied the radiation in participative medium. Nowadays, as the rapidly increase of computer power, the use of MCM in radiation calculation becomes more and more popular. Comprehensive reviews about the implementation of MCM in radiative transfer are available in refs. [START_REF] Howell | The Monte Carlo method in radiative heat transfer[END_REF][START_REF] Walters | Monte carlo methods for radiative heat transfer in scattering media[END_REF].

In the conventional MCM, a large number of photon bundles carrying a fixed amount of radiative energy are emitted in the system and their history is traced until the energy carried is absorbed at a certain point in the participative medium or at the wall, or until it exits the system. The departure point, propagation direction, spectral frequency of the ray and the point where the energy is absorbed are independently and randomly chosen according to given distribution functions. A main advantage of MCM is that it can be easily employed in complicated systems and can incorporate all the important physical effects (such as the anisotropic scattering distribution, spectral dependence of wall or medium properties, the directional dependence of wall properties) without any assumption and a large increment of computional effort. Moreover, the statistical features of the results of MCM enables the system uncertainty, commonly represented by the standard deviation, to be computed in addition to other quantities. This is generally unavailable in other more deterministic approaches. However, the drawback attendant is the need of a large number of rays to obtain statistically meaningful results. In fact, the standard deviation tends to be proportional to 1/ √ N [START_REF] Howell | The Monte Carlo method in radiative heat transfer[END_REF], where N is the total number of bundles. Therefore, an increment of three times of ray number only reduces the uncertainty by one half.

Different approaches for variance reduction

In order to reduce the long convergence time and the large memory requirement of MCM, several approaches, such as the "energy partitioning" method, the reverse MCM and the reciprocal MCM, have been proposed in the literature.

• "Energy partitioning" method

In conventional MCM, the energy carried by each ray is fixed until either the energy is totally absorbed at a single point or it exits the domain. This method, referred to as "standard ray tracing" method here, is very inefficient in an open configuration with highly reflective walls and/or optically thin medium, since most of the rays exit the system without any contribution to the statistical but consuming CPU time [START_REF] Modest | Determination of radiative exchange factors for three dimensional geometries with nonideal surface properties[END_REF][START_REF] Shamsundar | Monte carlo solutionseffect of energy partitioning and number of rays[END_REF]. To alleviate this obstacle, the "energy partitioning" method [START_REF] Walters | Monte carlo methods for radiative heat transfer in scattering media[END_REF] is proposed, where the energy carried by a ray is attenuated gradually by absorption of the cells along the optical path until its depletion or it leaves the enclosure. Each ray in this method can provide more statistical samples than the one in the "standard ray tracing" method and a smaller standard deviation is resulted. Consequently, this method is more efficient, particularly in the open system aforementioned, a tremendous computer time savings is achieved [START_REF] Modest | radiative heat transfer[END_REF]. Due to its high efficiency, this "energy partitioning" method is applied in this thesis.

• Reverse MCM

The solution of the entire enclosure is calculated simultaneously in conventional MCM (also called forward MCM) by emitting in all the cells. Generally, this method is reasonably efficient for a majority of systems where the overall radiation field is desired. However, if only the radiative intensity hitting a small detector and/or over a small solid angle is needed, this method can become inefficient since the ray which hits the small detector is only a very small fraction of the emitted bundles inside a large enclosure while all the other rays have no contribution. In order to obtain an accurate solution of the intensity incident at the detector, a prohibitive large number of photon bundles are needed, resulting in a huge computational cost. To alleviate this problem, the reverse MCM (also called backward MCM) based on the principle of reciprocity, can be applied instead. As shown in Fig. 1.4 (b), in order to calculate P ea AB (the radiative power emitted by radiation source A and absorbed by the detector B), the bundle is followed in a reverse direction, from the termination site (the detector B) to the participative medium A, rather than track it from A to B as in forward MCM ( Fig. 1.4 (a)). Hence, in reverse MCM, only the rays having a contribution to the solutions have been followed, and this greatly enhances the efficiency [START_REF] Walters | Rigorous development for radiation heat transfer in nonhomogeneous absorbing, emitting and scattering media[END_REF][START_REF] Walters | Monte carlo methods for radiative heat transfer in scattering media[END_REF]. [START_REF] Li | Importance of turbulence radiation interactions in turbulent diffusion jet flames[END_REF] has used the reverse MCM to compute the apparent directional emissivity of one-dimensional absorbing-emitting-scattering semitransparent slab and it is revealed that the reverse MCM converges more quickly than the forward MCM. However, generally, a large radiation source is required in this reverse MCM.

In order to extend its use into a system with collimated irradiation or point sources, Modest (2003a) has improved it by separating direct and scattered radiation.

• Reciprocal MCM Another two well known difficulties of conventional MCM are: 1) the treatment of quasi-isothermal system is difficult as the energy balance is very small compared to emission and absorption. A small error in the prediction of emitted and/or absorbed energy may cause a catastrophic error of the energy balance; 2) the computational effort of conventional MCM rapidly increases with the optical thickness [START_REF] Howell | The Monte Carlo method in radiative heat transfer[END_REF]. In fact, if the medium is optically thick, the mean optically free path is short, most of the rays emitted are directly absorbed in the vicinity of the cell under consideration and only a very few bundles participate in the distant Note that both of these two problems are associated with the violation of reciprocity principle. A main interest of a reciprocal MCM is to alleviate these problems since this method intrinsically fulfills the reciprocity principle. As shown in Fig. 1.4 (c), instead of using the path between A and B only for calculation of P ea AB (energy emitted by A and absorbed by B), the same path is used for the exchanged energy which is equal to the difference between P ea AB and P ea BA (energy emitted by B and absorbed by A). Here P ea AB is calculated directly while P ea BA is computed by using the reciprocity principle. In reciprocal MCM, the energy transport from A to B and from B to A are considered simultaneously. This reciprocal MCM was applied by [START_REF] Cherkaoui | Monte Carlo simulation of radiation in gases with a narrow-band model and a Net-Exchange Formulation[END_REF][START_REF] Cherkaoui | Radiative net exchange formulation within one-dimensional gas enclosures with reflective surfaces[END_REF] to simulate radiation in the limited cases of an one-dimensional quasi-isothermal slab with diffusely or specularly reflecting surface and non-gray gases. They reported that this method converged at least two orders faster than the conventional Monte-Carlo method and was much less sensitive to optical thickness. Moreover, in order to insure satisfactory convergence of reciprocal MCM for all types of optical thicknesses but always in 1D slabs, De [START_REF] De Lataillade | A net-exchange Monte Carlo approach to radiation in optically thick systems[END_REF] applied an adapted optical path sampling procedure which makes the Reciprocal MCM remain efficient even for very high optical thicknesses. Tésse et al. have also compared two different approaches of the Reciprocal MCM with the conventional MC method in one-dimensional benchmark cases with participative gases [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF]. These methods have been applied to 3D cases in Ref. [START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF]). The comprehensive formulations of the reciprocal approaches will be detailed later in section 1.2.2.

In addition to these approaches, biasing techniques have also been implemented. For example, Ref. [START_REF] Martin | Monte carlo analysis of the backscattering of radiation from a sphere to a plane[END_REF] have successfully used the biasing of sampled directions towards the parts of the system that most contribute to the addressed radiative quantity, while the biasing of sampled frequencies as function of temperature field and spectral properties are used in Ref. [START_REF] Dufresne | Inverse gaussian k-distributions[END_REF]).

Formulations of Reciprocal Monte-Carlo method

Any complex system, which consists of non-isothermal walls and non-isothermal, heterogeneous, emitting-absorbing participative medium, can be discretized into N v and N s isothermal, homogenous finite cells of volume V i and surface S i respectively. The radiative power P i (W ) in any cell i can be written as the (1.20)

For volume cells, for instance, P exch ij is defined as:

P exch ij = P ea ji -P ea ij = ∞ 0 (P ea ν,ji -P ea ν,ij )dν, (1.21)
where P ea ν,ij is the spectral radiative power emitted by volume i and absorbed by volume j. The exchanged power is calculated by using the reciprocity principle in reciprocal MCM.

Reciprocity principle

As shown in Fig. 1.5, dP ea ν,ij , the energy emitted by the differential volume element dV i and absorbed by dV j (=dA j × ds j ), is the energy emitted from dV i , transmitted by the media and absorbed within dA j ds j , or:

dP ea ν,ij = [4π κ ν (T i ) I 0 ν (T i ) dV i ] × dA j 4πr 2 × (τ ν,r ) × [κ ν (T j ) ds j ] (1.22)
where κ ν the spectral absorption coefficient. I 0 ν (T ) is the equilibrium spectral intensity. r is the distance between dV i and dV j and τ ν,r is the corresponding spectral transmissivity. This equation can be recast as

dP ea ν,ij I 0 ν (T i ) = τ ν,r κ ν (T i ) κ ν (T j ) dV i dV j r 2 (1.23)
Similarly, we can also demonstrate that dP ea ji can be expressed as

dP ea ν,ji I 0 ν (T j ) = τ ν,r κ ν (T i ) κ ν (T j ) dV i dV j r 2 (1.24)
Note that the right sides of these two Eqs. (1.23) and (1.24) are the same, so we obtain

dP ea ν,ij I 0 ν (T i ) = dP ea ν,ji I 0 ν (T j ) , (1.25)
and this is the so-called "reciprocity principle" which indicates that the ratio between dP ea ν,ij and dP ea ν,ji is equal to the corresponding equilibrium spectral intensity ratio. The "reciprocity principle" can be also applied to energy transfer between surface elements or between surface and volume elements. It is also valid in scattering or reflecting media.

Formulation of exchanged power

Using the reciprocity principle Eq. (1.25), the corresponding exchanged power between dV i and dV j can be expressed as [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF]:

dP exch ν,ij = dP ea ν,ji -dP ea ν,ij = τ ν,r κ ν (T i ) κ ν (T j ) [I 0 ν (T j ) -I 0 ν (T i )]
dV i dV j r 2 . (1.26) If we substitute dV j with the solid angle dΩ i and ds j by using the following relationship

dV j r 2 = ds j dΩ i (1.27)
then Eq. (1.26) yields

dP exch ν,ij = τ ν,r κ ν (T i ) κ ν (T j ) ds j [I 0 ν (T j ) -I 0 ν (T i )] dV i dΩ i (1.28)
An integration of dP exch ν,ij gives

P exch ij = +∞ 0 κ ν (T i )[I 0 ν (T j ) -I 0 ν (T i )] V i 4π A ij ν dΩ i dV i dν, (1.29)
where (1.30) and N p represents the total number of crossing of the cell j by a given optical path issued from the cell i. For the m-th crossing of the given optical path, τ ν (BF m ) is the spectral transmissivity between the source point B in cell i and the inlet point F m in the cell j (as shown in Fig. 1.6). α j m is the spectral absorptivity associated with the column m, defined as

A ij ν = Np m=1 τ ν (BF m ) α j m ,
α j m = 1 -exp [-κ ν (T j )l jm ] (1.31)
where l jm is the length of m-th column. The exchanged radiative energy between cell i and cell j can be calculated from Eq. (1.29) where the reciprocity principle is clearly fulfilled.

"thesis_yufang" -2013/11/27 -19:31 -page 38 -#38

✐ ✐ ✐ ✐ ✐ ✐ 38 1.2. MONTE-CARLO METHOD B dV i V j V i dV j F 1 F 2 l j1 l j2

Probability Density Functions

In MCM, a lot of photon bundles are emitted in the enclosure in order to calculate the radiative quantities. The departure point B, propagation direction ∆ and wave-number ν of these bundles are determined randomly according to the corresponding Probability Density Function (PDF) which will be presented in this part.

In fact, by introducing the emitted power P e i (T i ) of V i , which is defined as

P e i (T i ) = 4π V i +∞ 0 κ ν (T i )I 0 ν (T i )dν, (1.32)
Eq. (1.29) can then be transformed into

P exch ij = P e i (T i ) +∞ 0 [ I 0 ν (T j ) I 0 ν (T i ) -1] V i 4π A ij ν f i (B, ∆, ν)dΩ i dV i dν, (1.33)
where f i (B, ∆, ν) is a possible joint PDF, expressed as

f i (B, ∆, ν) dV i dΩ i dν = κ iν I 0 ν (T i ) dV i dΩ i dν P e i (T i )
.

(1.34)

Since the position B, direction ∆ and wave-number ν are statistically independent, f i (B, ∆, ν) is separated into three parts: f V i , f ∆i and f νi , given

f i (B, ∆, ν) dV i dΩ i dν = f V i (B)dV i f ∆i (∆)dΩ i f νi (ν)dν = 1 V i dV i 1 4π dΩ i κ ν (T i )I 0 ν (T i ) +∞ 0 κ ν (T i )I 0 ν (T i )dν dν (1.

39

Here the three independent function f V i (B), f ∆i (∆) and f νi (ν) can be used as the PDF for B, ∆ and ν respectively, given

f V i (B) = 1 V i , (1.36) f ∆i (∆) = 1 4π , (1.37) f νi (ν) = κ ν (T i )I 0 ν (T i ) +∞ 0 κ ν (T i )I 0 ν (T i )dν
.

(1.38) Furthermore, since the direction ∆ is determined by the azimuthal angle φ and the polar angle θ, the solid angle Ω i is replaced with θ i and φ i , and we obtain

f ∆i (∆)dΩ i = 1 4π i dΩ i = sin(θ i ) 2 dθ i 1 2π dφ i = f θi (θ i )dθ i f φi (φ i )dφ i . (1.39)
Here the PDF for φ and θ are respectively given as:

f θi (ν) = sin(θ i ) 2 , (1.40) f φi (ν) = 1 2π . (1.41)
All of these PDFs are independent and they will be applied for the determination of photon bundles.

Determination of a bundle

In order to mimic the radiation correctly in MCM, the bundles should be statistically emitted according to the PDFs. Assuming that N i bundles carrying equal energy P e i /N i are emitted from cell V i , the probability that a given bundle has a wave-number between ν and ν + dν is equal to f νi (ν)dν. Moreover, the probability that a given bundle has a wave-number between 0 and ν is given by the cumulative distribution function R ν , defined as

R ν = ν 0 f νi (ν)dν = ν 0 κ ν (T i )I 0 ν (T i )dν +∞ 0 κ ν (T i )I 0 ν (T i )dν
.

(1.42)

It indicates that if we want to simulate the emitted energy from V i by using N i bundles of equal energy, the fraction R ν (ν) of these bundles should have a wavenumber lower than ν. To realize this, a random number uniform distributed between 0 and 1 can be used. In fact, the probability that a uniform distributed random number has a value R ν0 is also R ν0 . For a given random number R ν0 , the corresponding wave-number ν 0 can be computed implicitly from Eq. (1.42 number R by using Eq. (1.42), then the fraction of the wave-number smaller than ν 0 is R ν0 and the cumulative distribution function is fulfilled. Similarly, to choose the polar angle θ and azimuthal angle φ, a random number R θ and R φ are generated respectively, given

R θ = θ 0 f θ (θ i )dθ i = θ 0 sin(θ i ) 2 dθ i , (1.43) R φ = φ 0 f φ (φ i )dφ i = φ 0 1 2π dφ i . (1.44)
Then the angles can be explicitly calculated from these two independent uniform distributed random number R θ and R φ as

θ = arccos(1 -2R θ ), (1.45) φ = 2πR φ . (1.46)
For an opaque surface node with diffuse emissivity, the distribution function for θ changes to

R θw = θw 0 f θw (θ i )dθ i = θw 0 sin(θ i )cos(θ i )dθ i , (1.47)
and angle θ w is calculated as

θ w = arccos( R θw ).
(1.48)

Generally, as discussed in Ref. [START_REF] Zhang | Radiation Monte Carlo Approach dedicated to the coupling with LES reactive simulation[END_REF], the departure point B can also be randomly chosen by following three steps: i) create a smallest parallelepiped surrounding cell V i ; ii) randomly pick one point within this parallelepiped; iii) if the picked point is within cell V i , then use this point as the departure point B. If not, repeat the second step until the picked point is inside cell V i . However, in Emission-based Reciprocity Method, an "infinitesimal volume" dV i centered in a given point B can be applied and now all the photon bundles are directly started from B, avoiding the process to randomly choose the departure point. This will be discussed later in detail.

Statistical estimations

After being emitted from the source point, each photon bundle crosses successively each cell in the selected optical direction. When the bundle goes through a cell, the energy absorbed in this cell is computed with the local absorptivity (based on the crossed length) and the remained energy in the path leaving from this cell can be computed with the local transmissivity. In this way, the energy carried by the bundle will decrease gradually until the amount of energy becomes less than a cutoff value or until the bundle leaves the enclosure. Finally, for a given shot n, the spectral radiative power emitted by a cell i that is absorbed by a cell j writes [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF])

P ea ijn ν = P e i (T i ) N i A ijn ν , (1.49)
where N i is the total number of bundles from cell i. According to the reciprocity principle, the radiative power which is exchanged between i and j, for the shot n, is

P exch ijn ν = P ea jin ν -P ea ijn ν = P ea ijn ν [ I 0 ν (T j ) I 0 ν (T i ) -1].
(1.50)

Note that here P ea ijn ν is calculated in a deterministic way (use Eq. (1.49) while P ea jin ν is computed from the reciprocity principle as

P ea jin ν = P ea ijn ν I 0 ν (T j ) I 0 ν (T i )
.

(1.51)

If a large number of photon bundles N i are emitted from cell i, after tracing the history of these bundles, the statistical estimation of P exch ij and P ea ij can be calculated by averaging the contribution of all the optical paths as:

P ea ij = N ij n=1 P ea ijn ν
(1.52)

P exch ij = N ij n=1 P exch ijn ν (1.53)
where N ij is the number of optical paths which give rise to a non-zero contribution and • denotes a statistical estimation value.

The conventional forward MCM is based on P ea ij . In addition, by using P exch ij , two reciprocal MCMs, the Emission-based Reciprocity Method (ERM) and Absorption-based Reciprocity Method (ARM), are also proposed in Ref. [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF]. The principles of these three models are detailed in the following part.

Forward Method (FM)

As mentioned previously, in FM, a photon emitted from cell V j and arriving at cell V i is only used for the energy transport from V j to V i (P ea ν,ji ). By summing the contribution of all the rays started from all the cells V j of the enclosure (as shown in Fig. 1.7), the radiative power of cell V i is statistically estimated as

P F M i = Nv+Ns j=1 P ea ji -P e i (T i ), (1.54)
where P e i (T i ) is calculated deterministically from Eq. (1.32). Obviously, the calculation of P F M i depends on the emission of all the cells in the system.

Emission-based Reciprocity Method (ERM)

In the Emission Reciprocity Method (ERM), the statistic estimation of the radiative power in the cell i is calculated as

P ERM i = Nv +Ns j=1 P exch ij .
(1.55)

As shown in Fig. 1.7, in ERM, all the photon bundles are started from cell i, which indicates only the bundles from the cell V i are needed for the computation of P ERM i . Therefore, for the situation where only the solution of a certain domain is needed, this method is efficient. Moreover, by controlling the bundle number emitted from each cells, the uncertainty can be locally controlled. This feature is useful if highly accurate solution is only required in a certain part of the system. Note that, for cell i, the emitted energy is calculated in a deterministic way while the absorbed energy is computed by using the reciprocity principle. Obviously, the accuracy of deterministic calculation is higher than reciprocal calculation. Therefore, the computed emitted energy is more accurate than absorbed energy and hence ERM is more adapted to the zone where emission is dominant than absorption (high temperature zone).

Absorption-based Reciprocity Method (ARM)

In ARM, a lot of bundles are shot from every cell of the system. By following the progress of all the bundles and using the reciprocity principle, the radiative power in cell V i is calculated as

P ARM i = -

43

more bundles are needed in every cell of the entire domain. By contrast to ERM, for cell V i , the absorbed energy is deterministically calculated whereas the reciprocity principle is applied for emitted energy. Consequently, ARM is more adapted for the low temperature region.

In addition to these two reciprocal methods, an Optimized Reciprocity Method (ORM) is proposed by [START_REF] Dupoirieux | An optimized reciprocity monte carlo method for the calculation of radiative transfer in media of various optical thicknesses[END_REF] where ERM is used in the case when T i is higher than T j , and ARM in the opposite case.

"Infinitesimal" emission source volume

In ERM, by substituting Eq. (1.50) and (1.49) into Eq. (1.53), then the radiative power of cell V i can be expressed as

P ERM i = P e i (T i ) N i Nv+Ns j=1 A ijn ν [ I 0 ν (T j ) I 0 ν (T i ) -1] = P e i (T i ) S.
(1.57) where

S = 1 N i Nv+Ns j=1 A ijn ν [ I 0 ν (T j ) I 0 ν (T i ) -1].
(1.58)

And the mean volume radiative power P ERM i (W/m 3 ) of cell V i can be calculated as

P ERM i = P ERM i V i = P e i (T i ) V i S = 4π S +∞ 0 κ ν (T i )I 0 ν (T i )dν (1.59)
where V i disappears on the right side. It implies that, the mean volume radiative power P ERM i is not affected by the volume of the emission source in ERM. Hence, instead of V i , we can apply an "infinitesimal volume" dV i centered at point B to calculate P ERM i (B) at position B, and now all the optical paths can be directly emitted from B rather than randomly choose the departure point inside V i . Moreover, the isothermal and homogenous assumption in the "infinitesimal volume" is now physically true. If assume that the volume power is a constant inside V i , the mean volume power of V i is then equal to P ERM i (B). However, this "infinitesimal" emission source volume is only valid for ERM. In fact, for ARM as an instance, the radiative power of V i is calculated as

P ARM i = Nv+Ns j=1 P e j (T j ) N j A jin ν [ I 0 ν (T i ) I 0 ν (T j ) -1] = Nv+Ns j=1 4πV j +∞ 0 κ ν (T j )I 0 ν (T j )dν N j A jin ν [ I 0 ν (T i ) I 0 ν (T j ) -1].
(1.60) and the volume power is

P ARM i = P ARM i V i = 1 V i Nv+Ns j=1 4πV j +∞ 0 κ ν (T j )I 0 ν (T j )dν N j A jin ν [ I 0 ν (T i ) I 0 ν (T j ) -1] (1.61)
Obviously, P ARM i depends on the emission source volume V j and the "infinitesimal" volume can not be used.

Statistical Variance

Since Monte Carlo methods are statistical methods, the solution will generally fluctuates randomly around the correct answer and these fluctuations will decrease as the number of samples increases. In order to estimate the error associated with the statistical result X(N ), the N photon are broken up into I subsamples and normally, each subsample would include identical amounts of bundles, leading to

N i = N I i = 1, 2, . . . , I.
(1.62) 63) where X(N i ) is the mean value of the i-th subsample. The I subsamples may then be treated as if they were independent experimental measurements of the same quantity. And the corrected sample standard deviation can be calculated as

X(N ) = 1 N N i=1 N i X(N i ) = 1 I n i=1 X(N i ) i = 1, 2, . . . , I. (1.
σ 2 = 1 (I -1) I i=1 X(N i ) -X(N ) 2 i = 1, 2, . . . , I.
(1.64)

while the variance or corrected standard deviation of the sample mean X(N ) equals

σ 2 m = σ 2 I .
(1.65)

Based on the central limit theorem, we can say with 68.3% confidence that the correct answer lies within the limits of X(N ) ± σ m , with 95.5% confidence within X(N ) ±2σ m , or with 99% confidence within X(N ) ± 2.58σ m . The determination of the variance associated with a computed quantity is an important advantage of the MCM, when compared to other more deterministic models for which the accuracy is not easy to determined.

Introduction

Monte Carlo method is widely used in many applications involving radiative transfer [START_REF] Howell | The Monte Carlo method in radiative heat transfer[END_REF][START_REF] Farmer | Comparison of monte carlo strategies for radiative transfer in participating media[END_REF]. At the beginning it was mainly used to produce reference solutions to validate more approximated techniques (see for instance, [START_REF] Guo | Radiative heat transfer in inhomogeneous, nongray, and anisotropically scattering media[END_REF][START_REF] Marakis | A parametric study of radiative heat transfer in pulverised coal furnaces[END_REF][START_REF] Mishra | Computational efficiency improvements of the radiative transfer problems with or without conduction -a comparison of the collapsed dimension method and the discrete transfer method[END_REF]). But the method is now directly used for determining the effects of radiation in combustion media [START_REF] Zhang | Monte Carlo References method of radiative transfer applied to a turbulent flame modeling with LES[END_REF][START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF][START_REF] Lataillade | Modelisation detaillee des transferts radiatifs et couplage avec la cinetique chimique dans des systemes en combustion[END_REF], for studying turbulence-radiation interactions [START_REF] Coelho | Numerical simulation of the interaction between turbulence and radiation in reactive flows[END_REF][START_REF] Wu | Direct numerical simulation of turbulence/radiation interaction in premixed combustion systems[END_REF][START_REF] Snegirev | Statistical modeling of thermal radiation transfer in buoyant turbulent diffusion flames[END_REF]), scattering media [START_REF] Stankevich | Monte Carlo ray-tracing simulation of light scattering in particulate media with optically contrast structure[END_REF] or collimated pulse irradiation [START_REF] Wu | Monte Carlo simulation of transient radiative transfer in a medium with a variable refractive index[END_REF], etc. Walters et al. [START_REF] Walters | Rigorous development for radiation heat transfer in nonhomogeneous absorbing, emitting and scattering media[END_REF][START_REF] Walters | Monte carlo methods for radiative heat transfer in scattering media[END_REF] have proposed a reverse Monte Carlo algorithm based on the reciprocity principle. This approach generally accelerates the convergence of the Monte Carlo method. Cherkaoui et al. [START_REF] Cherkaoui | Monte Carlo simulation of radiation in gases with a narrow-band model and a Net-Exchange Formulation[END_REF][START_REF] Cherkaoui | Radiative net exchange formulation within one-dimensional gas enclosures with reflective surfaces[END_REF] (2002), allow the determination of the radiative power in any cell for these complex systems. These authors have shown that ERM and ARM are more suitable than the Forward Monte Carlo Method (FM) for optically thick or nearly isothermal media. ARM is efficient in the calculation of radiative power or flux in low temperature regions, while ERM is, on the contrary, more efficient in hot regions. [START_REF] Dupoirieux | An optimized reciprocity monte carlo method for the calculation of radiative transfer in media of various optical thicknesses[END_REF] have proposed an optimized reciprocity method (ORM) by combining ERM and ARM. However, to compute the radiative power in some precise cells, FM, ARM or ORM requires to generate optical paths from all the cells of the system, while ERM is only based on the optical paths issued from the cells where results are expected. Therefore, ERM is more efficient when the calculation is limited to a part of the system. Consequently the aim of the present paper is to improve the efficiency of ERM, in particular in the treatment of the radiative fluxes and powers in the colder regions. The resulting Optimized Emission Reciprocity Method (OERM) is based on a frequency distribution function associated with the system maximum temperature. Section 2.3 briefly summarizes different approaches of the reciprocity method.

The principles of the ERM optimization are defined in Sec. 2.4 and first applied to a simple case involving a grey medium. Results of OERM for four benchmark cases, typical of combustion gases, are discussed in Sec.2.5. Comparisons with ERM results of Ref. [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF] are also carried out.

Different reciprocity Monte Carlo approaches

Any complex system can be discretized into N v and N s isothermal finite cells of volume V i or surface S i , respectively. The radiative power in any cell i can be written as the sum of the exchange powers P exch ij between i and all the other cells j, i.e. (2.1)

For volume cells, for instance, P exch ij is given by

P exch ij = +∞ 0 κ ν (T i ) [I • ν (T j ) -I • ν (T i )] V i 4π A ij ν dΩdV i dν, (2.2)
where I • ν (T ) is the equilibrium spectral intensity and κ ν (T i ) the spectral absorption coefficient relative to the cell i. A ij ν accounts for all the paths between emission from any point of the cell i and absorption in any point of the cell j, after transmission, scattering and possible wall reflections along the paths. Its expression is detailed in Ref. [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF] given for exchanges between a volume cell and a surface or between two surface cells. In a Monte Carlo method, a large number of shots N i , each of them characterized by a source point, a direction and a frequency ν, are stochastically generated from a cell i. For a given shot n, the spectral radiative power emitted by a cell i absorbed by a cell j writes

P ea ijn ν = P e i (T i ) N i A ijn ν , (2.3)
where P e i (T i ), the total emitted power of cell i and it is expressed as:

P e i (T i ) = 4πV i +∞ 0 κ ν (T i )I • ν (T i )dν.
(2.4)

According to the reciprocity principle, the radiative power which is exchanged between i and j, for the shot n, is

P exch ijn ν = P ea ijn ν I 0 ν (T j ) I 0 ν (T i ) -1 (2.5)
Statistical estimations P ea ij and P exch ij of P ea ij and P exch ij are obtained, at the large number limit, by summing the contributions of all the shots which connect i and j, i.e.

P ea ij = N ij n=1
P ea ijn ν , (2.6)

P exch ij = N ij n=1 P exch ijn ν , (2.7)
where N ij is the number of the shots which are issued from cell i and have a contribution to P ea ij or P exch ij . Accurate statistical determination of P exch ij can only be obtained if the number of shots connecting the cells i and j is large. For complex 3D geometrical configurations, characterized for instance by about 10 4 cells, a huge 10 4 × 10 4 matrix has then to be generated for each frequency. The previous approach is only realistic for simple configurations, typically for one-dimensional systems. Consequently, simplified reciprocity Monte Carlo approaches have been developed for complex systems [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF][START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF][START_REF] Dupoirieux | An optimized reciprocity monte carlo method for the calculation of radiative transfer in media of various optical thicknesses[END_REF]) by cumulating during the computation the results related to a cell i. In the Emission-based Reciprocity Method(ERM) and Absorption-based Reciprocity Method (ARM), the radiative power in a cell i are respectively estimated by (2.9) without storage of P exch ij ν . The principles of ERM and ARM are summarized in Fig. 2.1. P ERM i and P ARM i can be simultaneously computed. An optimization can be realized [START_REF] Dupoirieux | An optimized reciprocity monte carlo method for the calculation of radiative transfer in media of various optical thicknesses[END_REF]) by using ERM when the exchange between cell i and cell j is controlled by the error of cell i meaning that T i > T j , and ARM in the opposite case. But ARM, as the Forward Monte Carlo Method (FM) which is classically defined as

P ERM i = Nv+Ns j=1 P exch ij , (2.8)
P F M i = Nv+Ns j=1 P ea ji -P e i (T i ), (2.10)
presents an important drawback compared to ERM. In these two approaches, paths have to be shot from all the cells of the system, even if the results are only required in some critical cells. On the contrary, by using ERM, paths are only shot from cells where the results are wished. Moreover, when radiation is coupled to turbulent transfer, generally in reacting systems, much more accurate computations are required close to an emission point than far from this point. The use of ERM is also in this case more pertinent than using ARM.

Consequently, an Optimized Emission-based Reciprocity Method (OERM) is developed in the next section.

Optimization principle

2.4.1 Analysis of ERM limitations

P exch ij (Eq. (2.
2)) can be expressed as (2.11) where the distribution functions of frequency, volume and solid angle are [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF])

P exch ij = P e i (T i ) +∞ 0 I • ν (T j ) I • ν (T i ) -1 × V i 4π A ijν f ν i (ν)dν f V i dV i f Ω i dΩ i ,
f ν i (ν) = κ ν (T i )I • ν (T i ) +∞ 0 κ ν (T i )I • ν (T i )dν ; f V i = 1 V i ; f Ω i = 1 4π .
(2.12)

At the limit of large shot number, statistical estimation P exch ij is as In a Monte Carlo method, the frequency ν n associated with a shot n is generally obtained by using a uniform random number R n belonging to the range [0, 1]; ν n is then determined from the implicit equation

P exch ij = P e i (T i ) N i N ij n=1 I • νn (T j ) I • νn (T i ) -1 A ijn νn (2.
R n = νn 0 f ν i (ν)dν = νn 0 κ ν (T i )I • ν (T i )dν +∞ 0 κ ν (T i )I • ν (T i )dν . (2.14)
The limitations of this approach applied to ERM are illustrated by a simple example for which an analytical reference solution [START_REF] Taine | Transferts thermiques : Introduction aux transferts d'énergie[END_REF] (in Fig 2 .4). The frequency distribution function used in ERM is based on the spectral emitted power. Consequently, as shown in Fig. 2.3, most of the shots issued from cells at 700 K are characterized by low frequencies. But the power that is absorbed by a cold cell at 700 K has mainly be emitted by hot regions, characterized by much higher frequencies. P a ν is then strongly underestimated by ERM in the case of cells at 700 K. This phenomenon obviously does not appear for hot cells at 2,500 K as the emitted radiation spectrum is very close to the absorbed one, as shown in Fig. 2.4. Moreover if emission at low frequencies by cells at 700 K is underrepresented in the spectrum at 2,500 K, it weakly contributes to the absorbed power of a cell at 2,500 K.

OERM

The key idea of this optimized method is to use a frequency distribution function independent of the local emission distribution function. This frequency distribution function has to be suitable for both emission, that characterizes hot regions, and absorption. A simple choice is the emission distribution function which is associated with T max , the highest temperature within the system, i.e:

f ν (ν, T max ) = κ ν (T max )I • ν (T max ) +∞ 0 κ ν (T max )I • ν (T max )dν
.

(2.15)

In these conditions, the radiative exchange power between i and j given by Eq.

(2.2) can be expressed as

P exch ij =P e i (T max ) +∞ 0 I • ν (T i ) I • ν (T max ) κ ν (T i ) κ ν (T max ) I • ν (T j ) I • ν (T i ) -1 f ν (ν, T max )dν f V i dV i f Ω i dΩ i , (2.16)
The optimized Monte Carlo method, based on N i shots n of frequencies ν n from the cell i, allows now P exch ij to be estimated as

P exch ij = P e i (T max ) N i N ij n=1 I • νn (T i ) I • νn (T max ) κ νn (T i ) κ νn (T max ) I • νn (T j ) I • νn (T i ) -1 A ijn νn
(2.17) and finally P OERM i is computed by using Eq. (2.8). It is the Optimized Emission-based Reciprocity Method (OERM), characterized by the frequency distribution function defined by Eq. (2.15). Note that, if the frequency distribution function is associated with the emission distribution function at the maximum temperature, a corrective factor Tmax) is applied to each shot in order to rigorously obtain the emission distribution function at temperature T i in a cell i. The results of OERM applied to the example of Sec.2.4.1, in the same computation conditions, are shown in Fig. 2.2. They agree very well with the reference solution for the whole temperature range, contrary to the case of ERM results, according with the previous analysis of Sec.2.4.1. Indeed, the standard deviation is nearly uniform along y, while in the ERM case, it is much higher in the cold zone near the wall. Note that, in the hot central region of the slab, the results of ERM and OERM are very close: as the two frequency distribution functions are practically identical, OERM turns into ERM in this region. The averaged wall flux obtained from OERM is 556.4 kW/m 2 with a standard deviation of only 16.4 kW/m 2 . This result is very close to the reference value of 559.6 kW/m 2 , contrary to the result of ERM (418.3 kW/m 2 ). [START_REF] Goody | The correlated-k method for radiation calculations in nonhomogeneous atmospheres[END_REF], by using the parameters of [START_REF] Soufiani | High temperature gas radiative property parameters of statistical narrow-band model for H 2 O, CO 2 and CO, and correlated-K model for H 2 O and CO 2[END_REF]. Cases 1 and 2 correspond to a thin medium: the total emissivities of the whole gas medium, assumed to be isothermal at 2500K and 500K, are equal to 0.03 and 0.15, respectively. Cases 3 and 4 correspond to a globally thick medium. For cases 1 and 3, the gas mixture is hot and the walls are cold; the contrary occurs for cases 2 and 4. Results obtained from OERM and ERM are compared. They are also compared with FM and ARM results of Ref. [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF]), when available. There is no physical discrepancies between ERM results of this work and of Ref. [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF]. Figures 2.5 and 2.8 deal with radiative power per unit volume for cases 1 and 3 and Table 2 with the wall flux. As expected, the average values of the radiative power per unit volume and the radiative flux obtained from OERM are always close to the reference ones. Figure 2.6 shows that the associated standard deviations obtained by ERM in all cases is much higher than that of OERM in cold zones, near a wall or at a wall, but similar results are found in the core of the gaseous medium by the two methods. Moreover, if the average flux value predicted by ERM is acceptable in case 1 (thin medium) when compared to the reference one, in the thick medium case (case 3) it strongly differs from the reference solution and has an important standard deviation. In the last case the Monte Carlo computation is far from the convergence. In these common cases (cold walls and hot gas), the frequency distribution function of OERM allows again the absorption by the cold regions to be accurately computed, contrary to the case of ERM. Figure 2.9 (for case 3) and Table 2 (for cases 1 and 3) also compare the previously obtained standard deviations with those issued from ARM and FM, in Ref. [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF]). It appears that the standard deviations issued from ARM and FM are often close to those obtained by OERM, however, in some cases, they are smaller. For the Forward Method (FM) and the Absorption- based Reciprocity Method (ARM), the explanation is simple: it is always the frequency distribution function of the cell emitting the shots which is used for calculating absorption by a cold region. However FM and ARM present many drawbacks detailed in Ref. [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF]. Moreover, by using these approaches, computations have to be completely achieved in all the cells of the system, and not only in cells where a result is needed. Figures 2.7 and 2.10 present radiative power per unit volume and Table 2 the wall flux, for cases 2 and 4. In these much less common cases, conclusions are similar to the previous ones. OERM leads to much better results than ERM in the cold region, i.e. in the central region of the gas mixture. The two models agree in the hot regions, close to the walls. As the walls are hot, there are only small discrepancies on the wall flux results obtained by OERM and ERM.

I • νn (T i ) I • νn (Tmax) κν n (T i ) κν n (

Conclusion

A complete Monte Carlo reciprocity method, based on complete calculation of exchange powers between cells is not realistic for system involving absorbingemitting gases in complex geometrical configurations. An Emission-based Reciprocity Method (ERM) allows paths of the Monte Carlo method to be shot only from the points where results are wished, which is a determinant computational advantage, by comparison with Forward Monte Carlo Methods (FM) or Absorption-based Reciprocity Methods (ARM). The main drawback of the classical approaches of ERM is its inaccuracy in the treatment of absorption by cold regions. This disadvantage has been eliminated by using a frequency distribution function based on the emission distribution at the maximum temperature encountered in the system, in an Optimized Emission-based Reciprocity Method (OERM). The real emission distribution function of a cell at any temperature is then obtained by applying a corrective factor to each shot.

Models and numerical approaches for the flow

Firstly, the general transport equations of fluid flows and the three numerical approaches for turbulent flows, DNS, LES and RANS, are discussed. Then, wall functions associated with boundary layers are presented. Finally, the code YALES2 used for flow field simulations and its coupling with the Monte-Carlo method are detailed.

Transport equations of fluid flows

As proposed by Taylor and Kármán in 1937 [START_REF] Goldstein | Modern Developments in Fluid Dynamics[END_REF], turbulence can be defined as "an irregular motion which in general makes its appearance in fluids, gaseous or liquid, when they flow past solid surfaces or even when neighboring streams of the same fluid flow past or over one another". It is characterized by its irregularity in time and space and continuous spectra of length and time scales. In most of the practical engineering applications, flows are turbulent since the Reynolds number is large [START_REF] Wilcox | Turbulent modeling for CFD[END_REF]). This turbulent nature of flows introduces important difficulties in its understanding and description. However, even though turbulence is chaotic, it is deterministic and is described by transport equations. The continuity equation, momentum and energy balance equations and the ideal gas equation of state for three-dimensional and turbulent flows with variable medium properties can be written in the following form: where ρ, p and u i are the fluid mass density, pressure and velocity components, respectively. r is the ideal gas constant per mass unit. The enthalpy per mass unit h is expressed:

∂ρ ∂t + ∂ρu i ∂x i = 0, (3.1) ∂ρu i ∂t + ∂ρu i u j ∂x j = - ∂p ∂x i + ∂τ ij ∂x j + S i i ∈ [1, 3], (3.2) ∂ρh ∂t + ∂ρu j h ∂x j = ∂p ∂t + ∂pu j ∂x j + τ ij ∂u j ∂x i - ∂q cd j ∂x j + Q, ( 3
h = ∆h 0 + T T 0 c p (T ′ ) dT ′
, where c p is the mixture thermal capacity at constant pressure, T the temperature, T 0 a reference temperature and ∆h 0 the corresponding standard formation enthalpy. The viscous shear stress tensor τ ij writes

τ ij = µ ∂u i ∂x j + ∂u j ∂u i - 2 3 µ ∂u k ∂x k δ ij , (3.5)
where µ is the dynamic viscosity. Only homogeneous gas compositions will be considered, the total heat flux is then only composed of the conductive heat flux q cd i that is given by Fourier's law as

q cd i = -λ ∂T ∂x i , (3.6)
where λ is the thermal conductivity. S i and Q are the source terms in momentum and energy balance equation, respectively. When radiation is taken into account, Q is equal to radiative power per unit volume P R which is defined in section 1.1.

Direct Numerical Simulation

In Direct Numerical Simulation (DNS), the Navier-Stokes equations are solved without any approximate modeling (as shown in Fig. 3.1). All the spatial and temporal scales of turbulence should be resolved on the computational grid. Therefore, this leads to severe requirements on space and time resolution.

Resolution Requirement

Regarding space discretization, since the smallest turbulent length scales are of the same order of magnitude as the Kolmogorov scale which is defined as η = (ν 3 /ǫ) 1/4 (here ν = µ/ρ is the kinematic viscosity and ǫ is the dissipation rate of turbulent kinetic energy), the grid spacing of DNS should have the order of Kolmogorov scale O(η). Nevertheless, in order to ensure the relevant large scale to be represented properly, the computational domain should be large enough. In inhomogenous direction, the domain size is determined by the geometrical constraints while in homogeneous direction, the velocity fluctuation should be uncorrelated within half of the domain size.

For time advancement, explicit scheme or semi-implicit scheme are often applied. The time steps for explicit time advancement are generally smaller than the Kolmogorov time scale t η = (ν/ǫ) 1/2 from the linear stability criteria. In the near wall region, the viscous time scale ν/u 2 τ (u τ is the friction velocity on the wall) is also used to limit the time steps [START_REF] Friedrich | Direct numerical simulation of incompressible turbulent flows[END_REF]). In semiimplicit time advancement, the time step could be larger. Commonly implicit time advancement is only used for viscous term while explicit time advancement is still applied for convection term (semi-implicit time advancement) since the use of implicit time advancement for convection term is generally not suitable [START_REF] Moin | Direct Numerical Simulation: A Tool in Turbulence Research[END_REF].

Application of DNS

Although DNS has a huge requirement on resolution, it remains a valuable tool in turbulence studies. DNS results are used to analyze the turbulence physics, to assess turbulent modeling or even to quantify the accuracy of experimental techniques. DNS has been performed in various frameworks such as homogeneous turbulence [START_REF] Vincent | The satial structure and statistical properties of homogeneous turbulence[END_REF][START_REF] Jiménez | The References 219 structure of intense vorticity in isotropic turbulence[END_REF]), mixing-layers [START_REF] Vreman | Compressible mixed layer growth rate and turbulence characteristics[END_REF][START_REF] Maghrebi | DNS of forced mixing layer[END_REF], pipe flow [START_REF] Eggels | Fully developed turbulent pipe flow: a comparison between direct numerical simulation and experiment[END_REF][START_REF] Feldmann | Direct numerical simulation of fully developed turbulent and oscillatory pipe flows at Re τ =1440[END_REF], boundary layers [START_REF] Manhart | DNS of a turbulent boundary layer with separation[END_REF][START_REF] Na | Direct numerical simulation of a separated turbulent boundary layer[END_REF], backward-facing step [START_REF] Le | Direct numerical simulation of turbulent flow over a backward-facing step[END_REF][START_REF] Barri | DNS of backward-facing step flow with fully turbulent inflow[END_REF], channel flow (Kim et al. 1987;[START_REF] Antonia | Low-reynolds-number effects on near-wall turbulence[END_REF]) and even combustion systems [START_REF] Wu | Direct numerical simulation of turbulence/radiation interaction in premixed combustion systems[END_REF][START_REF] Luo | Direct Numerical Simulation of Pulverized Coal Combustion in a Hot Vitiated Co-flow[END_REF]). Among them, turbulent channel flow is extremely useful for the study of wall-bounded flow since its geometry is simple and it has the fundamental nature to investigate the turbulence structures in the near wall region.

DNS of channel flow

In the pioneering work of Kim and Moin (Kim et al. 1987;Kim and Moin 1987), they have investigated both the turbulent statistics of velocity field (Kim et al. 1987) and the transport of passive scalars (Kim and Moin 1987) in a fully developed channel flow with a Reynolds number of 3300. Many other studies have been carried out for both compressible and incompressible channel flows.

• Compressible channel flow [START_REF] Coleman | A numerical study of turbulent supersonic isothermal-wall channel flow[END_REF] have performed DNS of a supersonic channel flow with isothermal walls in order to investigate compressibility effects in boundary layers. They have shown that such compressibility effects are mainly caused by the mean density variation and the Morkovin's hypothesis is almost valid while the Van Driest transformation is successful. Morkovin's hypothesis is that: for compressible flow at moderate freestream Mach numbers (M<=5), dilatation is small and any deviation from incompressible turbulence can be accounted for by mean variations of fluid properties. This is the basis for the Van Driest transformation, a velocity scaling that accounts for the fluid-property variations to collapse compressible flow data onto the "universal" incompressible velocity profile. The success of van Driest transformation has also been demonstrated in study of Huang et al. (1995). However, it is presented in Ref. [START_REF] Morinishi | Direct numerical simulation of compressible turbulent channel flow between adiabatic and isothermal walls[END_REF]) that Morkovin's hypothesis does not explain well the near-wall asymptotic behavior of wall-normal turbulence intensity.

The effect of compressibility is also investigated in details in various Refs. [START_REF] Foysi | Compressibility effects and turbulence scalings in supersonic channel flow[END_REF][START_REF] Heinz | Turbulent supersonic channel flow: Direct numerical simulation and modeling[END_REF]).

• Incompressible channel flow

The low-Reynolds number effect in the near-wall region of incompressible channel flows has been studied by [START_REF] Antonia | Low-reynolds-number effects on near-wall turbulence[END_REF] and [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Re τ = 590[END_REF]. It has been indicated that a significant low-Reynolds number effect appears in the near-wall region of flows with small Reynolds numbers. [START_REF] Kasagi | Direct numerical simulation of passive scalar field in a turbulent channel flow[END_REF] have carried out a detailed analyzes of turbulent statistics and budget of temperature variance, its dissipation rate and turbulent heat flux in a channel flow with iso-flux walls. Great attention has also been given to the Reynolds number effect [START_REF] Abe | Direct numerical simulation of a fully developed turbulent channel flow with respect to the Reynolds number dependence[END_REF]Abe et al. 2004;Abe et al. 2004;[START_REF] Kozuka | DNS of turbulent heat transfer in a channel flow with a high spatial resolution[END_REF]) and Pr number effect [START_REF] Kawamura | DNS of turbulent heat transfer in channel flow with low to medium-high Prandtl number fluid[END_REF]Abe et al. 2004;[START_REF] Kozuka | DNS of turbulent heat transfer in a channel flow with a high spatial resolution[END_REF][START_REF] Na | Use of direct numerical simulation to study the effect of Prandtl number on temperature fields[END_REF][START_REF] Schwertfirm | DNS of passive scalar transport in turbulent channel flow at high schmidt numbers[END_REF]) on turbulent statistics. Moreover, Matsubara et al. have performed DNS for the investigation of the spanwise heat flux in a channel flow with spanwise temperature gradient [START_REF] Matsubara | Spanwise Heat Transport in Turbulent Channel Flow With Prandtl Numbers Rang-References 221 ing From 0[END_REF]) while some other effect in channel flow, such as variable gas properties [START_REF] Nicoud | Numerical study of a channel flow with variable properties[END_REF], buoyancy (Garcia-Villalba and del Alamo 2011), chemistry [START_REF] Cabrit | Direct simulations for wall modeling of multicomponent reacting compressible turbulent flows[END_REF] widely studied. Additionally, modified channel configurations are used to examine some other factors, such as transverse curvature [START_REF] Neves | Effects of convex transverse curvature on wall-bounded turbulence. part 1. the velocity and vorticity[END_REF], transpiration [START_REF] Sumitani | Direct numerical simulation of turbulent transport with uniform wall injection and suction[END_REF] and rotation (Kristoffersen and Andersson 1993).

As the development of high-performance computation technology, DNS is applied in more and more complex configurations. However, since the computational resource requirement increases rapidly with Reynolds number (degrees of freedom increase as ≈ Re 9/4 ), its application is still limited to flows at low or medium values. In simulations of flows characterized with a high Reynolds number, RANS models are widely used and some of them will be presented in the next section.

RANS simulation

In RANS (Reynolds Averaged Navier-Stokes equations) approach, only the mean flow fields are resolved while all the turbulence scales are modeled (as shown in Fig. 3.1). The balance equations in RANS simulation are obtained by averaging the instantaneous transfer equations. Time-averaged continuity equation and momentum equation write

∂ρ ∂t + ∂(ρ u i ) ∂x i = 0, (3.7) ∂(ρ u i ) ∂t + ∂(ρ u i u j ) ∂x j = - ∂ρ u ′′ i u ′′ j ∂x j - ∂p ∂x i + ∂τ ij ∂x j + S i , (3.8)
where over-line • and tilde • denote Reynolds and Favre averaged quantities respectively in RANS. The Reynolds averaged molecular viscous stress term τ ij is generally approximated as

τ ij = µ ∂ u i ∂x j + ∂ u j ∂u i - 2 3 µ ∂ u k ∂x k δ ij . (3.9)
And -ρ u ′′ i u ′′ j is the so-called Reynolds stress tensor needed to be modeled. Note that the Favre average is used here instead of Reynolds average in order to avoid the mass source term ρ ′ u ′ i and other terms due to Reynolds average [START_REF] Poinsot | Theoretical and numerical combustion[END_REF]. To model the Reynolds stress tensor and close the averaged momentum equation, there are two different kinds of models: eddy viscosity models and Reynolds stress models.

"thesis_yufang" -2013/11/27 -19:31 -page 66 -#66

✐ ✐ ✐ ✐ ✐ ✐ 66 3.3. RANS SIMULATION

Turbulent viscosity models

Based on the Boussinesq assumption [START_REF] Boussinesq | Théorie de l'Écoulement Tourbillant[END_REF][START_REF] Tennekes | A First course in turbulence[END_REF], the Reynolds stress tensor is linked to the mean velocity gradient via the turbulent viscosity, given

-ρ u ′′ i u ′′ j = µ t ∂ u i ∂x j + ∂ u j ∂u i - 2µ t 3 ∂ u k ∂x k δ ij + 2 3 ρk, (3.10)
where the turbulent dynamic viscosity is expressed as µ t = ρν t (here ν t is the turbulent kinematic viscosity) and the turbulent kinetic energy k is defined as

k = 1 2 3 k=1 u ′′ k u ′′ k . (3.11)
In order to model the turbulent viscosity, three different approaches are proposed in the literature: algebraic model, one-equation model and two-equation model.

• Algebraic model After the introduction of the mixing length concept, an algebraic expression has been proposed by Prandtl where the turbulent viscosity is related to the mixing length l m as

ν t ≈ l 2 m | S|. (3.12)
Here S is the mean stress tensor, defined as

S = 1 2 ∂ u i ∂x j + ∂ u j ∂x i . (3.13)
This model is the so-called mixing-length model. One problem associated with this model is the unknown mixing length which depends on the flow. Other algebraic models are the Baldwin-Lomax model [START_REF] Baldwin | Thin-layer approximation and algebraic model for separated turbulent flows[END_REF] and the Cebeci-Smith model [START_REF] Cebeci | Analysis of Turbulent Boundary References 215 Layers[END_REF], which separate the flow into inner and outer layers and apply different algebraic models for each layer. These models are more frequently used in practical applications

• One-equation model An alternative to the algebraic model is the one-equation model where the turbulent viscosity is expressed as a function of turbulent kinetic energy k and a transport equation is solved for k. In one-equation models, the turbulent viscosity writes where C µ is the model constant and l is the length scale to compute. When the length scale is defined as l = k 3/2 /ǫ, C µ is generally set to 0.09.

ν t = C µ l √ k, ( 3
The main disadvantage of this model, as for the algebraic model, is that it is not possible to find a general expression for the length scale.

• Two-equation model In addition to the balance equation of k, another transport equation for the turbulent dissipation rate ǫ can be solved, leading to the k -ǫ model. The turbulent viscosity is then expressed as

ν t = C µ k 2 ǫ . (3.15)
Generally, the standard k -ǫ model is only suitable for flows far way from walls since fully turbulent conditions are assumed. The flow in the near wall region, characterized by a local low turbulent Reynolds number, is either modeled by using wall functions or solved by applying a low-Reynolds number k -ǫ model [START_REF] Patel | Turbulence models for nearwall and low reynolds number flows -a review[END_REF] which uses some dumping functions to make the terms in the balance equations behave in a right way when approaching the wall.

Besides k-ǫ model, other two-equation models such as k-ω (Wilcox 1993) model and k -τ model [START_REF] Speziale | A critical evaluation of two-equation models for near wall turbulence[END_REF]) have also been proposed where the balance equation for specific dissipation ω or turbulent time scale τ is derived respectively. Two-equation models work reasonably well for a large number of engineering flows in term of mean quantities. However, due to the isotropic eddy viscosity assumption, they are unable to account for curvature effects or irrotational strains.

Reynolds stress models

Reynolds stress models (RSM) avoid the isotropic eddy-viscosity hypothesis and close the Reynolds-averaged Navier-Stokes equations by solving transport equations for the transports of the Reynolds stress tensor. These models are more realistic than the eddy viscosity models. However, they introduce six additional equations describing Reynolds stresses which makes it hard to implement and to stabilize numerically. Generally, algebraic Reynolds stress models [START_REF] Wallin | An explicit algebraic reynolds stress model for incompressible and compressible turbulent flows[END_REF] solve algebraic equations for the Reynolds stresses, whereas differential Reynolds stress models [START_REF] Chen | Fundamentals of turbulence modeling[END_REF][START_REF] Cécora | Differential reynolds stress modeling for aeronautics[END_REF] solve differential transport equations individually for each Reynolds stress component.

Since RSM models the Reynolds stress tensor anisotropy in a more rigorous manner than the eddy-viscosity models, it has greater potential to give accurate predictions for complex flows. However, the fidelity of RSM predictions is still limited by the closure assumptions employed to model various terms in the "thesis_yufang" -2013/11/27 -19:31 -page 68 -#68

✐ ✐ ✐ ✐ ✐ ✐ 68 3.4. LARGE EDDY SIMULATION
exact transport equations for the Reynolds stresses. More details about the RSM are available in Ref. [START_REF] Launder | Advanced turbulence models for industrial applications[END_REF].

An intermediate approach between DNS and RANS is the Large Eddy Simulation (LES) where only part of the turbulent scales are modeled. LES will be presented in the following section.

Large Eddy Simulation

In large eddy simulation, only the large dynamically important scales are resolved while the contribution of the small, nearly isotropic turbulent scales are modeled by using Sub-Grid Scale (SGS) models (as shown in Fig. 3.1). Applying a filter operator, the instantaneous mass, momentum and energy balance equations yield

∂ρ ∂t + ∂(ρ u i ) ∂x i = 0, (3.16) ∂(ρ u i ) ∂t + ∂(ρ u i u j ) ∂x j = - ∂τ SGS ij ∂x j - ∂p ∂x i + ∂τ ij ∂x j + S i , (3.17)
∂(ρ h) ∂t + ∂(ρ u j h) ∂x j = - ∂q SGS j ∂x j + ∂p ∂t + ∂pu j ∂x j + τ ij ∂u j ∂x i - ∂q cd j ∂x j + Q, (3.18)
where over-line • and tilde • denote, in the case of LES, filtered and massweighted filtered quantities, respectively.

Based on the SGS eddy-viscosity concept, the unresolved SGS stress tensor τ SGS ij = -ρ( u i u j -u i u j ) and the turbulent SGS heat flux q SGS i = -ρ( hu j -h u j ) are modeled as

τ SGS ij - 1 3 δ ij τ SGS kk = -2ρν SGS (S ij - 1 3 δ ij S kk ), (3.19) q SGS j = -λ SGS ∂T ∂x j , (3.20)
where S ij is the filtered shear stress tensor and the SGS conductivity is given by

λ SGS = ρ c p ν SGS Pr SGS . (3.21)
There exist various turbulent SGS models in the literature, for instance, the Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: 1. the basic experiment[END_REF], the dynamic Smagorinsky model [START_REF] Germano | Turbulence: the filtering approach[END_REF][START_REF] Lilly | A proposed modification of the germano subgrid-scale closure method[END_REF][START_REF] Meneveau | A lagrangian dynamic subgrid-scale model of turbulence[END_REF], similarity model [START_REF] Bardina | Improved subgridscale models for large-eddy simulation[END_REF], mixed model [START_REF] Vreman | An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications[END_REF]), WALE model [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF] and σ model [START_REF] Nicoud | Using singular values to build a subgrid-scale model for large eddy simulations[END_REF]). Among them, the σ model, based on the "thesis_yufang" -2013/11/27 -19:31 -page 69 -#69

✐ ✐ ✐ ✐ ✐ ✐
Chapter 3 -Models and numerical approaches for the flow 69 singular value decomposition of the velocity gradient tensor, ensures the cube behavior in the near wall region and predict the classic periodic channel flow (which is studied in this thesis) with a slightly higher accuracy than the dynamic Smagorinsky model [START_REF] Nicoud | Using singular values to build a subgrid-scale model for large eddy simulations[END_REF]. Moreover, it is easy to implement and has a low computational cost. Therefore, the σ model is used for modeling ν SGS in this thesis and the SGS Prandtl number Pr SGS is set to 0.9. For a comprehensive information on SGS models, the readers are referred to reviews or books such as Refs. [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows[END_REF][START_REF] Pope | Turbulent Flows[END_REF][START_REF] Meneveau | Scale-invariance and turbulence models for large-eddy simulation[END_REF].

Limitation of wall-resolved LES

In order to resolve properly the large vortical motions, the filter width should be smaller than the inertial length scale by a minimal factor, for instance, in a channel flow, the ratio between the filter width ∆ and the local integral dissipation length scale L = k 3/2 /ǫ (here k is the turbulent kinetic energy and [START_REF] Baggett | Resolution requirements in large-eddy simulations of shear flows[END_REF]). In the near wall region of a wall-bounded flow, since the integral length scale becomes of the same order of magnitude as viscous scales, the grid resolution is high there and the total computational cost scales as Re 2.4 (Piomelli 2008). Hence, if the near wall region is well resolved in LES (named wall-resolved LES), most of the resources are then used by the inner layer which is only 10% of the flow. For a Reynolds number ≈ 10 4 , about 50% of the resources are used for the inner layer [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF]. Obviously, the high cost of resolving the inner layer makes the wall-resolved LES inappropriate for industrial applications characterized by high Reynolds number. In order to overcome this obstacle of LES, there are two types of approaches, namely wall-modeled LES and hybrid LES.

ǫ is the turbulent dissipation rate) is ∆/L ≈ 1/10 (

Wall-modeled LES

In wall-modeled LES, instead of resolving the momentum transport in the inner layer, an approximated boundary conditions are applied, where the wall shear stress is related to the outer layer velocity. An equilibrium wall-stress model or a two-layers model could be used to obtain approximated boundary conditions.

Equilibrium model

In the equilibrium model, it is assumed that in the inner layer, a constant shear stress layer exists which implies that the effect of pressure and acceleration term are omitted. Meanwhile, if the first grid point y 0 is located out of the viscous sub-layer, the molecular viscous term could be neglected, and hence the log-law yields

u + || = u || u τ = 1 κ log y o u τ ν + B, (3.22)
where u || is the mean longitudinal velocity at y 0 and u τ is the friction stress. κ is the von Kármán's constant while B is the loglaw intercept. ν is the kinematic viscosity.

The equilibrium model was firstly proposed by [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large reynolds numbers[END_REF] in a turbulent channel flow with infinite Reynolds number. However, since the grid in the outer LES region was too coarse, the prediction in this case was poor when compared to experimental results. Later in 1975[START_REF] Schumann | Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli[END_REF] successfully applied it in a calculation of turbulent channel flow and annuli flow where the shear stress τ at walls are computed as

τ xy,w (x, z) = τ w u(x, y o , z) u(x, y o , z), τ yz,w (x, z) = ν w(x, y o , z) y 0 (3.23)
where x, y and z are the streamwise, wall normal and spanwise direction respectively.

• denotes averaging over an xz plane. τ w is balanced by the imposed pressure gradient or is calculated iteratively requiring that the mean velocity u(x, y o , z) satisfies the logarithmic law (Eq. (3.22)) at y 0 . The results of Schumann's studies were in good agreement with reference data. However, since the original equilibrium models are based on the constant shear stress layer approximation, its application is limited to simple flows. In order to extend its use to more complex configurations, some modifications are carried out by considering the effect of the inclination of the elongated structure in near wall region [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF]), the pressure gradient [START_REF] Wang | LES with wall models for trailing edge aeroacoustics[END_REF][START_REF] Duprat | A wall-layer model for large-eddy simulations of turbulent flows with/out pressure gradient[END_REF], buoyancy [START_REF] Moeng | A large-eddy-simulation model for the study of planetary boundary-layer turbulence[END_REF] or chemistry [START_REF] Cabrit | Direct simulations for wall modeling of multicomponent reacting compressible turbulent flows[END_REF].

Two-Layers Model

An alternative to the equilibrium model is the two-layer model. This model has firstly been proposed by [START_REF] Balaras | Two-layer approximate References boundary conditions for large-eddy simulations[END_REF] by resolving the two-dimensional turbulent boundary layer equations on a fine embedded grid, given

∂u i ∂t = - ∂u i u j ∂x j - 1 ρ ∂p ∂x i + ∂ ∂x j (ν + ν t ) ∂u i ∂x j for i = 1, 3 u 2 = - y 0 ∂u 1 ∂x 1 + ∂u 3 ∂x 3 dy, (3.24)
where ∂p ∂x i is the pressure gradient at the first point of the outer layer. i = 1, 2 and 3 represents the streamwise, wall-normal and spanwise direction respectively. A simple mixing-length model was used for the turbulent viscosity ν t

ν t = (κy) 2 D(y)|S|, (3.25)
where |S| is the magnitude of strain rate and the expression of D(y) used in Ref. [START_REF] Balaras | Two-layer approximate References boundary conditions for large-eddy simulations[END_REF] is Results of this study show that LES with this two-layers model gave a good prediction of the flow in a plane channel, square duct or a rotating channel flow. This model was also extensively applied and assessed in different configurations [START_REF] Cabot | large-eddy simulation with wall models[END_REF][START_REF] Cabot | Approximate wall boundary conditions in the large-eddy simulation of high reynolds number flow[END_REF][START_REF] Wang | Dynamic wall modeling for large-eddy simulation of complex turbulent flows[END_REF][START_REF] Kemenov | Explicit small-scale velocity simulation for high-re turbulent flows[END_REF]. However, it is outlined in Refs. [START_REF] Cabot | Approximate wall boundary conditions in the large-eddy simulation of high reynolds number flow[END_REF][START_REF] Nicoud | Large eddy simulation wall-modeling based on suboptimal control theory and linear stochastic estimation[END_REF]) that since the filter width of the first few points near the wall is larger than the local integral length scale, sub-grid errors can be dominant in the near wall region.

D(y) = 1 -exp[-(y + /A + ) 3 ] (3.
Hence the information at the first off-wall point, provided from LES to the wall model, always suffers from sub-grid errors and also numerical errors. In order to increase the accuracy of the information fed from LES, [START_REF] Kawai | Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy[END_REF] proposed an effective strategy where, instead of the first off-wall point, a point located further away from the wall is used as the outer boundary of the wall model. It has been demonstrated that the prediction is improved when using this approach. Another way to increase the accuracy of the wall-modeled LES results is stochastic forcing which will be detailed later. Most of the models described so far are dedicated to the velocity field. In order to extend them to the temperature field, a thermal wall model has been proposed by [START_REF] Benarafa | Temperature wall modelling for large-eddy simulation in a heated turbulent plane channel flow[END_REF], given the equation for energy transport as

∂T ∂t = - ∂T u j ∂x j + ∂ ∂x j [(α + α t ) ∂T ∂x j ] + Q (3.27)
where T is the mean temperature and α is the molecular diffusivity. The turbulent diffusivity α t is calculated from the turbulent viscosity ν t and a prescribed turbulent Prandtl number. Banarafa et al. have applied this thermal wall model in a heated turbulent channel flow where the computed wall normal heat flux has shown to be accurate. This thermal wall model has also been successfully employed and modified by applying a different formula for turbulent Prandtl number in Ref. [START_REF] Rani | Boundary-Layer Equation-Based Wall Model for Large-Eddy Simulation of Turbulent Flows with Wall Heat Transfer[END_REF].

Hybrid LES

Another popular approach is the hybrid LES, where the simulation model is switched from RANS in the inner layer to LES in the outer layer, by a modification of a length scale [START_REF] Spalart | Comments of feasibility of les for wings, and on a hybrid RANS/LES approach[END_REF][START_REF] Nikitin | An approach to wall modeling in large-eddy simulations[END_REF] or the use of a blending function in the turbulent transport model [START_REF] Hamba | A Hybrid RANS/LES Simulation of Turbulent Channel Flow[END_REF][START_REF] Shur | A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities[END_REF].

In this approach, near the interface of RANS and LES, the resolved length scale is close to those characterized in RANS region which is commonly larger than the length scale of the outer LES region. Hence, the generation of resolved small eddy scales are delayed and this results in the 'logarithmic law mismatch' (the logarithmic sub-layer is displace upwards) [START_REF] Baggett | On the feasibility of merging LES with RANS for the near-wall regions of attached turbulent flows[END_REF]. In order to remove the logarithmic law mismatch, [START_REF] Piomelli | The inner-outer layer interface in large-eddy simulations with wall-layer models[END_REF] generation of resolved eddy scales. This artificial forcing has improved the hybrid LES results and the logarithmic law mismatch has been reduced. It has also been found by [START_REF] Radhakrishnan | Reynolds-averaged and large-eddy simulations of turbulent nonequilibrium flows[END_REF] that the hybrid LES performs well without any forcing mechanism in flows where the mean flow has some destabilizing perturbation. A more detailed description of the wall models and hybrid LES could be found in the reviews [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF][START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF][START_REF] Cabot | Approximate wall boundary conditions in the large-eddy simulation of high reynolds number flow[END_REF]. However, for numerical simulation of turbulent flows, a main issue always exists in the treatment of the near wall region. In the following section, boundary layer theory is detailed.

Boundary layer theory

When a fluid flows past a surface, or an object moves through a fluid, there is a thin layer in the immediate vicinity of the surrounding surfaces where the effects of viscosity are significant. This thin layer is the so-called boundary layer where the velocity rapidly increases from zero on the surface to the free stream value away from the surface. If the Reynolds number is low, the boundary layer has a laminar structure whereas for a high Reynolds number as encountered in most applications, the boundary layer becomes turbulent. A detailed description of the flow within the boundary layer is very important for many problems, including the skin friction drag on an object, the heat and mass transfer that occurs in combustors. In this section, non-dimensional mean flow quantities of boundary layer are presented. It is worth noting that following the Morkovin's hypothesis, development of boundary layer theory is commonly done by Reynolds averaged variables instead of Favre average. Effects of density fluctuation are therefore usually neglected.

Mean flows in wall units

As demonstrated in appendix A, in thin equilibrium turbulent boundary layer, such as the inner layer in a channel flow, the momentum and energy balance equation can be simplified as

dτ tot dy = d dy µ ∂u ∂y -ρ u ′′ v ′′ = 0, (3.28) dq tot dy = - d dy λ ∂T ∂y -ρ h ′′ v ′′ = 0, (3.29)
where τ tot and q tot are the total shear stress and heat flux respectively while ρ u ′′ v ′′ and ρ h ′′ v ′′ are the turbulent shear stress and turbulent heat flux respectively. Based on some assumptions, the dynamic and thermal boundary layer can be represented by wall functions as described in the following part.

Momentum boundary layer in wall units

An integration of the momentum equation (Eq. (3.28)) from the wall gives rise to

τ w = µ du dy -ρ u ′′ v ′′ (3.30)
where τ w = µ w du dy y=0 is the wall shear stress. The turbulent shear stress -ρ u ′′ v ′′ can be approximated with the Boussinesq assumption [START_REF] Boussinesq | Théorie de l'Écoulement Tourbillant[END_REF][START_REF] Schmitt | Large-eddy simulation and experimental study of heat transfer, nitric oxide emissions and combustion instability in a swirled turbulent high-pressure burner[END_REF]), given

-ρ u ′′ v ′′ ≈ µ t du dy , (3.31)
where the turbulent dynamic viscosity µ t can be modeled by using the Prandtl mixing-length approximation

µ t = ρl 2 m du dy . (3.32)
The mixing length l m in the near wall region linearly increases with the wall distance y as l m = κy where κ is the Von Kármán constant and usually takes the value 0.41 [START_REF] Pope | Turbulent Flows[END_REF]. Based on this, we obtain the expression for µ t as µ t = ρ(κy) 2 du dy .

(3.33) And Eq. (D.6) changes to

τ w = µ du dy + µ t du dy = µ du dy + ρ(κy) 2 du dy 2 . (3.34)
Using the gas properties on the wall, a non-dimensional wall normal distance y + and non-dimensional velocity u + can be defined as

u τ = τ w ρ w , y + = ρ w u w y µ w , u + = u u τ , (3.35)
and hence equation (3.34) is expressed as • Viscous sub-layer In the immediate vicinity of the wall where 0 < y + < 5, the effect of the turbulent shear stress is negligible compared to that of the molecular shear stress, and we have

µ µ w du + dy + + ρ ρ w (κy + ) 2 du + dy + 2 = 1. ( 3
µ µ w du + dy + = 1. (3.37)
After an integration, it yields

u + = y + 0 µ w µ dy ′+ . (3.38)
If the variation of µ is small, then a linear wall function is obtained (3.39) which indicates that the u + increases linearly with non-dimensional wallnormal distance y + in the very near wall region.

u + = y + ,
• Logarithmic sub-layer On the other hand, in the region where y + > 30 , the effect of the turbulent shear stress is predominant while the molecular shear stress is negligible. Then Eq. ( 3

u + V D = u + 0 ρ ρ w 1 2 du ′ + = 1 κ In(y + ) + C, (3.42)
where C is the logarithmic law intercept. Similarly, by omitting the variation of gas properties, the classic log law is obtained:

u + = 1 κ ln(y + ) + C. (3.43)
The Von Kármán constant κ and constant C are generally set to 0.41 and 5.0 respectively [START_REF] Pope | Turbulent Flows[END_REF]. However, for flows with a small Reynolds number, κ and C are commonly evaluated as 0.4 and 5.5 respectively (Kim et al. 1987). The transition layer between the Viscous sub-layer and the logarithmic sub-layer is the buffer layer. It is the location of the main production of turbulence.

Thermal boundary layer in wall units

Similarly, for the thermal boundary layer, an integration of the energy balance equation (Eq. (3.29)) leads to

q w = -λ ∂T ∂y + ρ h ′′ v ′′ . (3.44)
Based on the classic Reynolds analogy [START_REF] Reynolds | On the extent and action of the heating surface of steam boilers[END_REF]) between momentum and heat flux, the turbulent fluxes of enthalpy ρ h ′′ v ′′ can be modeled as

ρ h ′′ v ′′ ≈ ρ c p v ′′ T ′′ = -λ t dT dy (3.45)
where λ t is the turbulent heat diffusion coefficient, defined as

λ t = µ t c p Pr t = ρ c p (κy) 2 Pr t du dy . (3.46)
The turbulent Prandtl number Pr t here is commonly set to a constant as 0.85-0.9 in boundary layers [START_REF] Wilcox | Turbulent modeling for CFD[END_REF][START_REF] Kader | Temperature and concentration profiles in fully turbulent boundary layers[END_REF]. Introducing µ t , Eq. (3.44) can be rewritten as

q w = -λ dT dy -λ t dT dy = -λ dT dy - ρ c p (κy) 2 Pr t du dy dT dy . (3.47)
The non-dimensional temperature T + , based on the wall gas properties, as Similarly, for the viscous sub-layer and logarithmic sub-layer, Eq. (3.49) can be simplified and integrated as:

T + = |T -T w | T τ , T τ = |q w | c p w ρ w u τ , ( 3
• Viscous sub-layer

T + = - y + 0 c p w ρ w ν w λ dy ′+ . (3.50)
If the gas properties variation is omitted, the linear wall function is obtained as

T + = Pr y + .
(3.51)

• logarithmic sub-layer

T + 0 ρ ρ w 1 2 c p c p w dT ′+ = - Pr t κ In(y + ) + C ′ . (3.52)
Omitting the change in gas properties, it is obtained that

T + = - Pr t κ In(y + ) + C ′ . (3.53)
Note that the constant C' depends on the Prandtl number Pr. Based on experimental results, [START_REF] Kader | Temperature and concentration profiles in fully turbulent boundary layers[END_REF] proposed a logarithmic law for thermal boundary layer as

T + = α In(y + ) + β(Pr), (3.54)
where α = 2.12, that corresponds to Pr t ≅ 0.85, and the constant β is given as

β(Pr) = (3.85Pr 1/3 -1.3) 2 + 2.12In Pr.
(3.55)

Mean flows in semi-local coordinate

When the gas properties significantly vary in the flow, Huang et al. (1995) recommended to use the semi-local coordinates where the local gas properties are considered. This will be described in the following part.

Momentum boundary layer

For a momentum boundary layer, if the non-dimensional velocity and wall distance are defined based on the gas properties as Then for the viscous and logarithmic sub-layer, the wall functions obtained are:

u * τ = τ w ρ(y) , y * = ρ(y)u * τ y µ(y) , u * = u u * τ , ( 3
• Viscous sub-layer

u * = y * ; (3.58)
• logarithmic sub-layer

u * = 1 κ In(y * ) + C. (3.59)
Note that now the wall functions are "universal" for all thin equilibrium turbulent boundary inner layers.

Thermal boundary layer

Similarly, for thermal boundary layer, if define • logarithmic sub-layer

T * τ = |q w | c p (y) ρ(y) u * τ , T * = |T -T w | T * τ , ( 3
T * ≈ Pr t κ In(y * ) + C. (3.63)
It has been demonstrated by [START_REF] Dailey | Large eddy simulation of constant heat flux turbulent channel flow with property variations: Quasideveloped model and mean flow results[END_REF] that using semi-local coordinates, the velocity profiles nearly collapse for different cases of channel flow with variable density. The semi-local coordinates have also been supported by [START_REF] Coleman | A numerical study of turbulent supersonic isothermal-wall channel flow[END_REF] in their study of compressible channel flows.

The code YALES2

The unstructured parallel finite-volume code YALES has been developed by Dr. Vincent Moureau (CORIA, Rouen) since 2007. This is a versatile numerical solver for a broad range of flow problems (combustion, magnetohydrodynamics, multi-phase flows, etc.). It is capable to efficiently handle unstructured meshes of several billions of elements (Moureau et al. 2011a;[START_REF] Moureau | DNS analysis of a Re = 40,000 swirl burner[END_REF], thus enabling the DNS or highly resolved LES of laboratory and semi-industrial configurations. YALES2 is based on the Low Mach number approximation and the corresponding low-Mach Navier-Stokes equations are resolved with a projection method [START_REF] Chorin | Numerical solution of the navier-stokes equations[END_REF]) for constant or variable density flows.

Low-Mach Navier-Stokes equations

As demonstrated in Appendix B, based on the Low Mach number approximation implied in YALES2, the original balance equations of mass (Eq. (3.1)), momentum (Eq. (3.2)) and energy (Eq. (3.3)) and the ideal gas equation of state (Eq. (3.4)) can be written in the following form

∂ρ ∂t + ∂ρu i ∂x i = 0 (3.64) ∂ρu i ∂t + ∂ρu i u j ∂x j = - ∂p 2 ∂x i + ∂τ ij ∂x j + S i i ∈ [1, 3] (3.65) ∂ρh ∂t + ∂ρu j h ∂x j = ∂p 0 ∂t - ∂q j ∂x j + Q (3.66) p 0 ρ = rT (3.67)
where the pressure p is decomposed into two parts: a constant thermodynamic pressure p 0 and a pressure fluctuation p 2 . Based on this decomposition, the density is now directly calculated from p 0 in Eq. (3.67) rather than coupled with the p 2 .

Interest of Low Mach number approximation

In order to ensure the numerical stability when solving the governing equations, the time step needs to be limited by two restrictions namely, the Courant-Fredrichs-Lewy (CFL) condition and the restriction on the basis of grid-Fourier (Fo) numbers (for momentum, energy and scalars). Regarding the CFL condition, two number CFL a and CFL c , associated with the propagation of acoustic wave and convection respectively, are defined as where v is the velocity along the direction considered, c is the sound speed, and ∆x is the space step. When the diffusion terms in the governing equations are more important, the condition necessary to ensure stability is dictated by the restriction on the Fo numbers, for instance, the Fourier number for the energy equation is defined as

CFL a = (|v| + c)∆t ∆x , ( 3
Fo = α∆t ∆x 2 , (3.70)
where α is the thermal diffusivity. These characteristic numbers are used to limit the time step, where the smallest time step is chosen from these numbers.

In low Mach number, the limit of time step from the CFL a is the most serious [START_REF] Julien | Simulation numérique d'écoulements réactifs instationnaires à faibles nombres de Mach[END_REF], since

CFL a ∆x (|v| + c) < Fo ∆x 2 α < CFL c ∆x (|v|) . (3.71)
where the latter part of the above expression is generally valid for wall-resolved flows.

However, if the Low Mach number approximation is used, the propagation of acoustic wave disappears and hence the time step is released from the restriction of CFL a .

Numerics

For space discretization, a vertex-centered method is adopted in YALES2, where a dual set of control volumes is created, which are centered around the vertices of elements. For a detailed description of the spatial discretization in 3D, readers are referred to the Ph.D thesis of S. Vantieghem [START_REF] Vantieghem | Numerical simulations of quasi-static magnetohydrodynamics using an unstructured finite volume solver: development and applications[END_REF].

The spatial discretization available in YALES2 includes second order and fourth order centered finite-volume schemes. A detailed description of the higher order spatial discretization in available in Ref. [START_REF] Kraushaar | Application of the compressible and low-mach number approaches to large-eddy simulation of turbulent flows in aero-engines[END_REF]. To advance the solution in time, YALES2 disposes of explicit as well as implicit time integration schemes. The numerical schemes for time advancement in YALES2 are: 2ndorder or 4th-order Runge-Kutta scheme, Crank-Nicholson scheme and TFV4A scheme which is a blending of RK4 and TTG4A [START_REF] Moureau | Yales2 home page on www[END_REF]. In this thesis, the 4th-order centered schemes and 4th-order TFV4A explicit scheme are retained for spacial discretization and time advancement respectively.

Coupling between YALES2 and the Monte-Carlo method

In order to investigate the interaction of radiation and turbulent flows, it is necessary to fully couple the flow field code YALES2 and the Monte-Carlo method described in section 1.2.2. Several variables such as radiative power and temperature are exchanged on the interface of the two solvers while the synchronization in CPU time is also very important for the computational efficiency.

• Data exchange The procedure used for coupling YALES2 and the Monte-Carlo method is shown in Fig. 3.3. Firstly, to calculate the gas radiative properties and then the radiative field, the instantaneous temperature T and gas composition profile (mole fraction of species X i ) calculated in YALES2 are required in the Monte-Carlo method. However, in this thesis, since all the cases have a constant gas composition, only temperature profile is sent from YALES2 to the Monte-Carlo code. On the other hand, the radiative power term calculated in the Monte-Carlo code is required by YALES2. It acts as a source term in the balance equation of energy and hence modify the temperature field.

It is worth noting that the grid used in Monte-Carlo is not necessarily the same as the one for flow field. Hence, before sending the radiative power term or after receiving the temperature field, an interpolation is carried out in order to obtain the needed variable field over the grid.

• Synchronization in CPU time

In order to use the processors efficiently, one should pay attention to the synchronization in CPU time.

As discussed in Refs. [START_REF] Amaya | Unsteady coupled convection, conduction and radiation simulations on parallel architectures for combustion applications[END_REF][START_REF] Dos Santos | Coupled large eddy simulations of turbulent combustion and radiative heat transfer[END_REF], there is no need to couple the flow and radiation calculations every iteration. This is important, especially when Monte-Carlo method is applied, since the computation of radiation is much more expensive than the flow field calculation. In fact, as shown in Fig. 3.4, for a periodic channel flow which has a bulk Reynolds number of 12000 and two isothermal walls with wall temperature T w1 = 950K and T w2 = 1150K respectively (case C4 used in chapter 4), the change of instantaneous temperature profile during five iterations are very small. Therefore, it is reasonable to use the same radiative power field for five consecutive iterations of flow field. This indicates that for a coupled simulation, the radiation calculation is only performed every N i (N i = 5 for the case aforementioned) iterations of the flow field. This makes the computation much less expensive while ensuring the computation accuracy. Moreover, to avoid waiting of processors, one should make sure that: t r = N i • t f , where t r and t f are the time needed for one iteration of radiation and of flow field respectively. They depend on the number of processors used for each solver, i.e.

t f = γ f (P f ) t 1 f P f , (3.72) t r = γ r (P r ) t 1 r P r , (3.73)
where the subscript f and r denote the flow and radiation calculation respectively while γ(P ) is the speeding-up function. t 1 is the time required for one iteration on one processor and P is the number of processors. If the total number of processors used for the two solvers is P tot , a balanced distribution of the processors for flow field calculation is then

P f = P tot γ f t 1 f γ r t 1 r /N i + γ f t 1 f .
(3.74)

Conclusion

To simulate turbulent flows, there exist three different approaches: DNS, RANS and LES. DNS is the most accurate approach but its application is limited to flows at low or moderate Reynolds number due to its high resolution requirement. RANS models are commonly employed in engineering applications thanks to its low computational cost. However, the accuracy of its results are questionable in complex configurations. A technique intermediate between DNS and RANS is LES which models only the small turbulent scales. In LES of wall-bounded flows, the near wall region is generally modeled in order to alleviate the high cost to well resolve the boundary layer.

The boundary layer can be divided into two main layers and wall functions are retained for each layer. For flows with gas properties variation, the semi-local coordinate enables the obtained wall functions become "universal" by applying the local gas properties. A low-Mach finite-volume solver YALES2 is used for the flow field simulation of this work. And it is fully coupled with Monte-Carlo method to investigate radiation-flow field interaction. Special attention has been given to the data exchange of these two solvers and the synchronization in CPU time.

DNS of turbulent channel flow: effect of radiation on the mean flow fields

In this chapter, Direct numerical simulation of a CO 2 -H 2 O-N 2 mixture turbulent flow in a two-dimensional channel is performed. The fluid solver is coupled to a reciprocity Monte Carlo method based solver when radiation is taken into account. Effect of gas-gas radiation and gas-wall radiation on the mean temperature and flux field have been studied separately. The dependence of radiation effects on the set of conditions (temperature level, wall emissivity, pressure, Reynolds number) is also analyzed. A publication about this work is inserted: Y.F. Zhang, R. Vicquelin, O. Gicquel and J. [START_REF] Zhang | Optimized emission-based reciprocity monte carlo method to speed up computation in complex systems[END_REF]. Physical study of radiation effects on the boundary layer structure in a turbulent channel flow. International Journal of Heat and Mass Transfer 61 (2013) 654-666.

Additional information are available in annexes:

• The adequacy of computational domain of channel flow is validated in Appendix C.

• An original approach to specify the source term S i in the momentum equation in order to control the channel flow Reynolds number is detailed in Appendix D.

"

Abstract

A complete numerical coupling between radiation and turbulent convection in a channel gas flow has been performed for different temperature, optical thickness (pressure) and wall emissivity conditions. In this model, radiation is treated from the CK approach and a Monte Carlo transfer method; The flow by a Direct Numerical Simulation. Both the effects of turbulence on radiation fields and of radiation on turbulent fields are accounted for. Gas-gas and gas-wall radiation interactions generate antagonist effects on the temperature and flux fields. The first one tends to increase wall conductive flux while the second one to decrease it. Consequently, the structure of the temperature field and the wall conductive flux often strongly differ from results without radiation. Classical wall log-laws for temperature are then strongly modified by the global radiation effects. Many conditions encountered in applications are discussed in the paper. The observed modifications depend on all the set of conditions (temperature level, wall emissivity, pressure, Reynolds number), i.e. on the relative magnitudes of radiation gas-gas and gas-wall phenomena and of global radiation flux and conductive flux without radiation. Keyword: DNS, turbulence, gas radiation, Monte-Carlo, channel flow, temperature, wall-law

Introduction

In many industrial systems such as those dealing with combustion, conductive heat fluxes and radiative energy fluxes at walls condition the design stage and the material choice. Predicting these different fluxes with numerical simulations is therefore a great challenge that has been investigated in many works.

From the fundamental understanding of boundary layers and radiation energy transfer, models have been proposed to calculate this quantities. Although improving the prediction of heat fluxes remains an active area in each field, the determination of radiative fluxes and conductive fluxes at walls is always done separately without considering any coupling between these two modes of energy transfer. In this paper, the coupling of radiation and turbulent forced convection is studied to investigate whether the wall heat flux and the boundary layer structure can be modified by radiative energy transfer. The importance of Turbulence-Radiation Interaction (TRI) has been highlighted in several works. Experimental studies [START_REF] Gore | Spectral and total radiation properties of turbulent carbon monoxide/air diffusion flames[END_REF][START_REF] Ji | Thermal radiation properties of turbulent lean premixed methane air flames[END_REF][START_REF] Zheng | Measurements and inverse calculations of spectral radiation intensities of a turbulent ethylene/air jet flame[END_REF], theoretical analysis [START_REF] Burns | Turbulence radiation interaction modeling in hydrocarbon pool fire simulations[END_REF][START_REF] Baum | Radiation and velocity fields induced by localized temperature fluctuations[END_REF][START_REF] Soufiani | Radiation effects on turbulent heat transfer in channel flows of infrared active gases[END_REF]) and numerical simulations [START_REF] Wu | Direct numerical simulation of turbulence/radiation interaction in premixed combustion systems[END_REF][START_REF] Deshmukh | Direct numerical simulation of turbulence-radiation interactions in homogeneous nonpremixed combustion systems[END_REF][START_REF] Deshmukh | Direct numerical simulation of turbulence-radiation interactions in a statistically onedimensional nonpremixed system[END_REF][START_REF] Dos Santos | Coupled large eddy simulations of turbulent combustion and radiative heat transfer[END_REF]) have been carried out to investigate the effect of TRI in different systems and a comprehensive review is available in [START_REF] Coelho | Numerical simulation of the interaction between turbulence and radiation in reactive flows[END_REF][START_REF] Coelho | Turbulence-Radiation Interaction: From Theory to Application in Numerical Simulations[END_REF]. TRI is a consequence of the highly non-linear coupling between the radiative intensity and the turbulent temperature and gas species composition fields. It consists of two parts, namely the influence of turbulence on radiation and vice-versa. Regarding the former, turbulence leads to an increase in the medium transmissivity [START_REF] Jeng | Radiative heat fluxes near turbulent buoyant methane diffusion flames[END_REF][START_REF] Gore | Spectral and total radiation properties of turbulent carbon monoxide/air diffusion flames[END_REF], the radiative power [START_REF] Coelho | Detailed numerical simulation of radiative transfer in a nonluminous turbulent jet diffusion flame[END_REF][START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF]) and radiative heat loss [START_REF] Li | Importance of turbulence radiation interactions in turbulent diffusion jet flames[END_REF][START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF]). In order to isolate and quantify individual contributions to TRI in a statistically onedimensional premixed combustion system, Wu et al. ( 2005) have performed a direct numerical simulation (DNS) coupled with a Monte-Carlo method. The results reveal that the temperature self-correlation contribution is only dominant in the case with smallest optical thickness in their study. Among the three correlation terms of the absorption coefficient, the ones with the Planck function and with the incident radiative intensity are not negligible, even in the most optically thin case, while for a case with intermediate value of optical thickness, the three correlations were all significant. Deshmukh et al. [START_REF] Deshmukh | Direct numerical simulation of turbulence-radiation interactions in homogeneous nonpremixed combustion systems[END_REF][START_REF] Deshmukh | Direct numerical simulation of turbulence-radiation interactions in a statistically onedimensional nonpremixed system[END_REF]) have also studied different contributions to TRI with DNS in a statistically homogeneous isotropic non-premixed combustion system and a one-dimensional turbulent non-premixed flame. Only the latter one was fully coupled with radiation.

Regarding the effect of radiation on the flow, radiation interacts and modifies the temperature field in non-reactive flows [START_REF] Gupta | Large-Eddy Simulation of Turbulence-Radiation Interactions in a Turbulent Planar Channel Flow[END_REF][START_REF] Soufiani | Radiation effects on turbulent heat transfer in channel flows of infrared active gases[END_REF][START_REF] Ghosh | Effects of radiative heat transfer on the structure of turbulent supersonic channel flow[END_REF]) and in reactive flows [START_REF] Dos Santos | Coupled large eddy simulations of turbulent combustion and radiative heat transfer[END_REF][START_REF] Damien | Analysis of the interaction between turbulent combustion and thermal radiation using unsteady coupled LES/DOM simulations[END_REF]. In combustion applications, the change of maximum temperature has a significant effect on NO x emission [START_REF] Barlow | Scalar profiles and no formation in laminar opposed-flow partially premixed methane/air flames[END_REF]. Influence of radiation is not restricted to the average temperature field: The intensity of the temperature fluctuations changes when radiation is taken into account [START_REF] Damien | Analysis of the interaction between turbulent combustion and thermal radiation using unsteady coupled LES/DOM simulations[END_REF][START_REF] Soufiani | Temperature turbulence spectrum for high-temperature radiating gases[END_REF] Previous studies on the impact of radiation in a channel flow configuration [START_REF] Amaya | Unsteady coupling of Navier-Stokes and radiative heat transfer solvers applied to an anisothermal multicomponent turbulent channel flow[END_REF][START_REF] Soufiani | Application of statistical narrow-band model to coupled radiation and convection at high temperature[END_REF][START_REF] Soufiani | Radiation effects on turbulent heat transfer in channel flows of infrared active gases[END_REF][START_REF] Gupta | Large-Eddy Simulation of Turbulence-Radiation Interactions in a Turbulent Planar Channel Flow[END_REF][START_REF] Ghosh | Effects of radiative heat transfer on the structure of turbulent supersonic channel flow[END_REF] do not expose a general trend or understanding of the influence of radiation in turbulent boundary layers. In LES and RANS, wall models are necessary to predict wall heat flux in turbulent boundary layers. Among these models, standard wall laws are still widely used although they are only valid for relatively simple flows with zero-pressure gradient and constant fluid properties. In real applications with more complex conditions, standard wall laws might have a great inaccuracy. In order to extend their usage to practical systems, other physical effects need to be considered such as compressibility [START_REF] Foysi | Compressibility effects and turbulence scalings in supersonic channel flow[END_REF]Huang et al. 1995;[START_REF] Coleman | A numerical study of turbulent supersonic isothermal-wall channel flow[END_REF], Prandtl number effect [START_REF] Kader | Temperature and concentration profiles in fully turbulent boundary layers[END_REF][START_REF] Kawamura | DNS of turbulent heat transfer in channel flow with low to medium-high Prandtl number fluid[END_REF], streamwise pressure gradient (Huang and Bradshaw 1995;[START_REF] Nickels | Inner scaling for wall-bounded flows subject to large pressure gradients[END_REF] or chemical reaction [START_REF] Cabrit | Direct simulations for wall modeling of multicomponent reacting compressible turbulent flows[END_REF]. The objective of this paper is to study the radiation effects on the temperature field in turbulent boundary layers by considering a turbulent channel flow. DNS is performed for the flow field in order to generate high-fidelity data.

For radiation simulation, a reciprocal Monte Carlo method is employed. This method is generally regarded as the most accurate one and is widely used in many applications involving radiative transfer [START_REF] Coelho | Numerical simulation of the interaction between turbulence and radiation in reactive flows[END_REF][START_REF] Wu | Direct numerical simulation of turbulence/radiation interaction in premixed combustion systems[END_REF][START_REF] Zhang | Monte Carlo References method of radiative transfer applied to a turbulent flame modeling with LES[END_REF][START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF]). The two solvers are fully coupled when radiation is taken into account. The gas radiative properties are determined by means of the correlated k-distribution (CK) model or its weak absorption limit [START_REF] Soufiani | High temperature gas radiative property parameters of statistical narrow-band model for H 2 O, CO 2 and CO, and correlated-K model for H 2 O and CO 2[END_REF], depending on the pressure condition. A detailed analysis of gas-gas and gas-wall radiation effect is carried out in this paper. For the gas-gas radiation, only the energy exchange between a gas cell and the surrounding gas is accounted for, while gas-wall radiation includes only radiative energy transfer between a gas cell and the walls. A priori these two parts have different effects on the gas temperature profile and, hence, the wall heat flux. The numerical models and approaches are presented in Section 4.3. In Section 4.4, the effects of gas-gas and gas-wall radiation are studied for both large and intermediate optical thickness medium. Results associated with different wall emissivity, wall temperature and Reynolds number are discussed in Section 4.5.

Models and numerical approaches

Flow simulation

Under the low-Mach number approximation, the mass, momentum, energy balance equations and the ideal gas equation of state write respectively, in tensorial notations

∂ρ ∂t + ∂(ρu i ) ∂x i = 0, ( 4
∂(ρu i ) ∂t + ∂(ρu i u j ) ∂x j = - ∂p ∂x i + ∂τ ij ∂x j + S i , (4.2) ∂(ρh) ∂t + ∂(ρu j h) ∂x j = ∂p ∂t - ∂q cd j ∂x j + P R (4.3) p = ρ r T, (4.4)
where ρ, p and u i are the fluid mass density, pressure and velocity components, respectively. The enthalpy per unit mass h is expressed:

h = ∆h 0 + T T 0 c p (T ′ ) dT ′
, where c p is the mixture thermal capacity at constant pressure, T the temperature, T 0 a reference temperature and ∆h 0 the corresponding standard formation enthalpy. The viscous shear stress tensor τ ij and the conductive flux vector q cd i write

τ ij = µ ∂u i ∂x j + ∂u j ∂u i - 2µ 3 ∂u k ∂x k δ ij , (4.5
)

q cd i = -λ ∂T ∂x i , (4.6)
where µ is the dynamic viscosity and λ the thermal conductivity. S i is a forcing source term required in the following channel flow computations. It is uniform and is not null only in the streamwise direction. It plays the same role as the pressure gradient to drive the flow against viscous forces. P R is the radiative power per unit volume.

A finite-volume method for massively parallel computations on complex grids, suitable for variable density flow (code YALES2 (Moureau et al. 2011b;Moureau et al. 2011c)) is used to solve Eqs. (4.1)-(4.4) with associated boundary conditions. The code is a low Mach-number solver: the pressure p is split into a uniform thermodynamic pressure p 0 and a hydrodynamic one p 1 . p 0 is used in the equation of state to compute the gas density ρ, while p 1 is obtained by solving a Poisson equation to enforce mass conservation. Spatial gradients are calculated with a centered fourth order accurate scheme. Advancement in time is done with a fourth-order Runge Kutta scheme for velocity. For scalars (here enthalpy), the latter temporal scheme is blended with another fourth-order two-step scheme based on Taylor expansion [START_REF] Kraushaar | Application of the compressible and low-mach number approaches to large-eddy simulation of turbulent flows in aero-engines[END_REF]) that is more dissipative. The blending factor is set as small as possible and yet high enough to avoid spurious oscillations in the scalar field with centered schemes. In this study, a typical value of 0.05 is retained.

Non-dimensional quantities in channel flows

For the channel flow defined in with

ρ b = 2δ 0 ρdY 2δ , u b = 2δ 0 ρudY 2δ 0 ρdY , µ b = µ(T b ) ; (4.8)
The bulk enthalpy per mass unit and the friction velocity, based on the wall viscous stress τ w , by

h b = 2δ 0 ρuhdY 2δ 0 ρudY , u τ = τ w ρ w 1/2 , (4.9)
where indices b and w are related to bulk and wall quantities, respectively, and • refers to average quantities. The bulk temperature T b is evaluated from the bulk enthalpy : h b = h(T b ). When the channel flow case is not symmetrical (T w,c = T w,h ), the friction velocity is different on both walls: u τ,c = u τ,h and, therefore: Re τ,c = Re τ,h . The non-dimensional distance y + , based on the distance to a wall y, the nondimensional streamwise velocity u + , and the non-dimensional temperature T + are written

y + = ρ w y u τ µ w , u + = u u τ , T + = |T -T w | T τ , (4.10)
where T τ is the friction temperature defined by

T τ = |q cd w | ρ w c pw u τ , (4.11)
where q cd w is the wall conductive flux, chosen positive in the Y direction. According to this definition and Fig. 4.1, q cd w is negative on both walls since T w,c ≤ T w,h .

Validation of the flow simulation

In order to test the suitability of the code for direct numerical simulations, the channel flow case of Kim and Moin (Kim et al. 1987;Kim and Moin 1987) is computed and taken as reference for cases where temperature acts as a passive scalar. This case, named C0, is defined in Tab. 4.1. The discretization of the domain, defined in Tab. 4.2, is uniform along X and Z. ∆X, the cell length in the X direction, is expressed in wall units:

∆X + = ρ w ∆Xu τ µ w .
(4.12)

∆Y + and ∆Z + , associated with the cell lengths ∆Y and ∆Z in the Y and Z directions respectively, are defined similarly. The mesh is refined along the Y direction close to the wall where ∆Y + = 0.8. The size of the domain was chosen large enough for results to be independent of this parameter.

For this case C0 only, density and gas properties are uniform and T w,c = T w,h . Three temperature fields associated with three values of Prandtl numbers Pr 1 = 0.1, Pr 2 = 0.71 and Pr 3 = 2.0 are simulated. In this case, the temperature behaves like a passive scalar and non-dimensional results do not depend on the actual value of T w,c . A uniform source term is added in the energy equation as in Ref. (Kim and Moin 1987) to counterbalance conductive heat fluxes at the walls.

The non-dimensional mean velocity u + and the scalar T + are plotted in Fig. 4.2 and compared to numerical results of Kim and Moin (Kim et al. 1987;Kim and Moin 1987), and to the log-law obtained by these authors: u + = 2.5 lny + + 5.5, in Fig. 4.2 a and to the log-law of Kasagi et al. [START_REF] Kasagi | Direct numerical simulation of passive scalar field in a turbulent channel flow[END_REF]): T + = 2.78 lny + + 2.09, in Fig. 4.2 b. Good agreement is obtained between our results and those of the previous authors, demonstrating the accuracy of the chosen numerical set up to perform DNS of channel flows.

Radiation simulation

The general organization of the radiation model, based on a reciprocal Monte Carlo approach, has been detailed by Tésse et al. [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF]. This model has been previously applied to combusting media involving radiation [START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF]) and improved in [START_REF] Zhang | Optimized emission-based reciprocity monte carlo method to speed up computation in complex systems[END_REF]. Only the principles of the method are briefly summarized here.

In the approach of Ref. [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF], the radiation computational domain is discretized into N v and N f isothermal finite cells of volume V i and faces of area S i , respectively. The radiative power in any cell i is written as the sum of the exchange powers P exch ij between i and all the other cells j, i.e.

P i = Nv+N f j=1 P exch ij = - Nv+N f j=1 P exch ji . (4.13)
For volume cells, for instance, P exch ij is given by

P exch ij = +∞ 0 κ ν (T i )[I • ν (T j ) -I • ν (T i )] V i 4π A ij ν dΩ i dV i dν, (4.14)
where I • ν (T ) is the equilibrium spectral intensity and κ ν (T i ) the spectral absorption coefficient relative to the cell i. dΩ is an elementary solid angle and A ij ν accounts for all the paths between emission from any point of the cell i and absorption in any point of the cell j, after transmission, scattering and possible wall reflections along the paths. Its expression and similar expressions for exchanges between a volume cell and a surface cell or between two surface cells are detailed in Ref. [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF]). In the reciprocity Monte Carlo method, a huge number of optical shots are issued from the cells. Statistical estimation P exch ij of P exch ij are obtained by summing the contributions of the N ij shots that connect i and j, i.e.

P exch ij = N ij n=1 P exch ijn ν . (4.15)
In order to increase the computational efficiency, emission is here only carried out from arbitrary small spheres around grid points instead of using finite grid cells that are neither isothermal nor homogeneous. In this condition, the small emission spheres are isothermal, homogeneous and optically thin. A consequence of this choice is that the Absorption-based Reciprocity Method (ARM), which requires emission from finite cell, is not suitable. Only Emissionbased Reciprocity Methods (ERM) can then be used. In order to overcome some drawbacks of ERM compared to ARM in cold regions of the medium, the Optimized Emission Based Reciprocity Method (OERM), proposed by Zhang et al. [START_REF] Zhang | Optimized emission-based reciprocity monte carlo method to speed up computation in complex systems[END_REF], is here used for the radiation frequency treatment.

An other advantage of ERM and OERM method is to allow the Monte Carlo convergence to be locally controlled. In all the present simulations a radiative power standard deviation of 3% of radiative power maximum value has been imposed at any grid point (as shown in Fig. 4.3 for instance).

93

for atmospheric pressure condition, and by using the weak absorption limit for high pressure cases [START_REF] Soufiani | High temperature gas radiative property parameters of statistical narrow-band model for H 2 O, CO 2 and CO, and correlated-K model for H 2 O and CO 2[END_REF]). In the directions X and Z, periodicity conditions are also used for radiation simulation, i.e., if a shot exits the domain, for instance, at the point (L X , Y, Z), it will enter at the point (0, Y, Z) with the same propagation direction. The grid of the radiation model is three times coarser in X direction and two times in Y and Z directions than the corresponding grid of the flow model.

In the chosen DNS conditions, all averaged fluxes are uniform in directions X and Z. The averaged energy balance equation then writes, from Eq. ( 4.3),

d dY [ρ v ′′ h ′′ (Y ) + q cd (Y ) + q R (Y )] = d dY q tot = 0, (4.16)
where ρ v ′′ h ′′ , q cd , and q R and q tot are the averaged turbulent convective heat flux, conductive flux, radiative energy flux and total flux, respectively. This total flux is also uniform along Y . Note that, the radiative energy flux can be split into two parts :

q R (Y ) = q R * (Y ) + q R ww , (4.17)
where q R ww is the flux exchanged between the walls through the whole gaseous medium, that does not participate to the fluid energy balance and is zero if the walls are at the same temperature, as encountered in many applications. In the studied configurations, it is always uniform. Then Eq. (4.16) writes

ρ v ′′ h ′′ (Y ) + q cd (Y ) + q R * (Y ) = q w , (4.18)
where q w = q tot is the total flux exchanged between a wall and the gaseous mixture without the wall-wall radiation contribution. At the wall Y = δ, for instance, q w writes

q w = q cd (Y = δ) + q R * (Y = δ) = q cd w + q R * w . (4.19)
The balance of these three terms in Eq. (4.18) is presented in Fig. 4.4 for the case C2 without radiation and the related case C2R3 which accounts for radiative energy transfer. This illustrates that radiation modifies the balance of the terms in Eq. (4.19) and hence, the wall conductive flux, given in Tab. 4.4. These complex coupling effects are progressively detailed in the next sections.

Effects of Gas-Gas and Gas-Wall radiative interactions

The effects of radiation are due to different coupled phenomena, in particular Gas-Gas (GG) and Gas-Wall (GW) interactions, but also possible multiple reflections. In this section, the effects of gas-gas and gas-wall interactions are separately considered. In gas-gas interaction, a given gas volume only exchanges radiative energy with all other gas volumes, while the walls only partially reflect radiation. Emission by the walls is not accounted for in the GG calculations.

On the contrary, in gas-wall interaction, a given gas volume only exchanges energy with walls, while the other gas volumes only transmit energy. Emission by the other gas volumes is not accounted for in GW calculations. A coupled computational case that is limited to gas-gas (respectively gas-wall) radiative interactions is called CnRm_GG (respectively CnRm_GW). The obtained results are then compared with the corresponding cases CnRm, which include all the radiation interactions, and Cn, which does not take into account radiative energy transfer.

Large optical thickness medium

In the case C1R1, defined in Tables 4.1 to 4.3 and characterized by a high pressure, the optical thickness of the gaseous mixture is large: The global medium Hotell's transmissivity at 1000 K is equal to 0.271. This allows the effects of gas-gas and gas-wall radiative interactions to be studied with a weak radiative coupling between the two walls. Due to the small wall temperature difference, the variation in mass density in cases C1 and C1R1 is small. The profiles of the mean non-dimensional streamwise velocity u + , not shown here, are not different from the one of case C0.

The profiles of the mean temperature T associated with the cases C1 and C1R1 are compared in Fig. 4.5 a. The corresponding T + profiles, for both the cold and the hot sides, are plotted in Fig. 4.5 b. In the case C1, without radiation, the distribution of T is practically antisymmetric and the T + profiles are identical for the two sides. These results agree well with the results of Kim and Moin (Kim and Moin 1987), obtained under the assumptions that the temperature is a passive scalar. When only gas-wall radiation is considered, in the C1R1_GW case, the temperature gradient is smaller in the vicinity of a wall than in C1 case: The associated wall conductive fluxes presented in Table 4.4 are two to three times smaller than in C1 case. Indeed, the wall tends to impose its temperature to the fluid. Consequently, the temperature variation is higher in the core of the flow than in case C1. The conductive flux variations are more important at the hot side than at the cold one, as gas-wall radiative interactions strongly increase with the temperature.

On the other hand, in the case C1R1_GG, the gas-gas radiative transfer is a supplementary transfer that homogenizes the temperature field within the gas by comparison with the case C1, without radiation, as shown in Figs. 4.5 a and b . Consequently the temperature gradients and the conductive fluxes strongly increase at the two walls, as shown in Table 4.4. When all the radiative effects are accounted for, in the case C1R1, the gas-gas Reτ,c Re τ,h q cd w,c Table 4.4: Wall fluxes (in W/m 2 ) q cd w and q R* w for different cases at cold and hot walls (see Eqs. (4.16) and (4.17)). For each case CnRm, the relative variation of q cd w compared to the case Cn without radiation is put between parentheses. Flux values are rounded, typical errors are within 2-3 %. and gas-wall interactions, that have opposite effects, are coupled. As shown in Fig. 4.3, in the vicinity of the wall, the amplitude of the radiative power associated with gas-gas interaction is much larger than the one associated with gas-wall interactions. Consequently, in the present conditions, the wall temperature gradients and the conductive fluxes increase, as in the gas-gas case, at the two walls by comparison with the case C1 without radiation, as shown in Tab. 4.4 and Fig. 4.5 a.

q cd w,h q R* w,c q R* w,h (cold) (hot) (cold) (hot) (cold) (hot)
An important result is that, in this case, no log-law can be clearly identified.

Intermediate optical thickness medium

All the previous study has been again carried out for the case C2R1, characterized by a fluid at atmospheric pressure instead of 40 atm, but with the same data as for C1R1 (Re b , δ, T 1 , T 2 , ε 1 , ε 2 ). In these new conditions, the optical thickness of the medium is much smaller than previously. The global Hotell's transmissivity at 1000 K of the channel is equal to 0.811. Consequently, interaction phenomena between the two walls now occur, due to multiple reflections.

Results are shown in Figs. 4.6 a, b and c. The same analysis as previously can be drawn, but the effects related to gas-gas and gas-wall interactions are smaller than at high-pressure. As shown in Tab. 4.4 , the conductive fluxes obtained by only accounting for gas-gas interaction (case C2R1_GG) for the two sides are larger than the conductive fluxes associated with the case C2, without radiation. On the contrary, the conductive fluxes obtained in case C2 and by only accounting for gas-wall interaction (case C2R1_GW) are close, especially at the cold side. Finally, the conductive flux associated with all radiation effects (C2R1) is larger than in the C2 case at the two walls as in the optically thick case. Moreover, the temperature fields obtained by accounting for all radiation effects differ from the fields computed without radiation. Once again the profiles associated with the case C2 agree with the temperature profile of (Kim and Moin 1987) while it is not any more valid when radiation is accounted for. As for the case C1R1, the usual thermal log-law is not valid in this case.

Influences of different parameters

The radiation intensity field and, consequently, the radiative power and the temperature fields within the gaseous medium are strongly influenced by the wall emissivities, the wall temperatures and the bulk Reynolds number of the flow. The roles of these three quantities are studied in this Section.

Influence of wall emissivity

Three other couples of wall emissivities are now considered, in the conditions C2 of the flow at atmospheric pressure: These radiative conditions R2, R3 and R4 are defined in Tab. 4.3. From a practical point of view, an emissivity of The antagonist gas-gas and gas-wall effects also appear on the average radiative power field shown in Fig. 4.8 a. On both walls, the radiative power magnitude increases when wall emissivities decrease. This global effect can be analyzed from gas-gas and gas-wall exact contributions to the total radiative power in cases C2Rm. These contributions are given in Figs. 4.8 b and c. When wall emissivity decreases, the reflected flux and possibly the number of reflections increase. Consequently, the gas-gas interaction effects increase as seen in Fig. 4.8 b . On the other hand, the flux exchanged between the gas and the walls and, hence, the gas-wall interaction effects decrease as seen in Fig. 4.8 c. The shape of the gas-wall contribution to the total radiative power is explained by splitting the gas-wall interaction into a gas-cold-wall and a gashot-wall interactions. These latter contributions are given in Fig. 4.8 d. Finally, the global radiative power is the sum of all contributions and follows the same trend as the gas-gas interaction which is significantly larger than the gas-wall interaction in these cases.

The radiative conditions R4 (see Tab. 4.9 a, by comparing the cases C2, C2R3 and C2R4. For case C2R4, the hot black wall strongly imposes its temperature to the close gaseous layer. In the other hand, the reflecting cold wall strongly increase gas-gas radiative interactions and homogenizes the fluid temperature. The cumulation of these effects leads to: i) A much smaller temperature gradient and conductive flux at the hot wall than in case C2R3 associated with an emissivity equal to 0.1 (see Tab. 4.4 ); ii) A weak averaged temperature variation in the core of the gaseous medium; iii) Consequently the temperature gradient at the cold wall and the associated conductive flux are much larger than in case C2R3, as shown in Tab. 4.4 . The radiative power fields of the two extreme cases C2R3 and C2R4 differ in the vicinity of both walls as shown in Fig. 4.9 b. It can be explained by the previously discussed cumulative effects on the temperature profile : On the hot side, the gas temperature is close to the hot wall temperature, which decreases the magnitude of the radiative power; On the cold side, both the black hot wall and the large hot region contribute to increase the radiative power. Figure 4.10 shows the radiative power field in the hot half part of the channel. It appears that the radiative power in case C2R3 is dominated by gas-gas interaction. On the contrary, in the case C2R4, gas-wall interaction overcomes gas-gas interaction except for the close vicinity of the hot wall. For intermediate emissivity value, intermediate radiative power fields are encountered. Note that, in the present simple case as in most of the previous ones, both gas-wall and gas-gas radiative interactions are strongly modified by the wall reflection, as the optical thickness of the mediums is weak. Gas-gas and gaswall interactions are not isolated phenomena, as in case C1Rm at high-pressure, characterized by a large global optical thickness. The wall reflection law and the medium optical thickness have important effects on the temperature fields as they can enhance or reduce the gas-gas or the gas-wall contributions.

Influence of the temperature

In order to study the influence of both the temperature level and the temperature difference, the C3 and C3R1 cases are considered. They are characterized by wall temperatures equal to 950K and 2050K. On the cold side, these data correspond to classical combustion applications at high pressure. The associated results are compared with those of C1 and C1R1 cases. Contrary to the previous cases, the gaseous medium density strongly varies for C3 and C3R1 cases and variations in the averaged velocity profiles are now observed in Fig. 4.11 a. Similarly, differences between the u + profiles are observed in Fig. 4.11 b, showing that the classical velocity log-law is not valid anymore. However, if the Van Driest transformation, defined by

u + V D = u + 0 ρ ρ w du + , (4.20)
is used to account for variable density effects, all velocity profiles collapse to the usual log-law (see Fig. 4.11 c).

The averaged temperature profile of C3R1 case is compared in Fig. 4.12 a with those of C3, C1 and C1R1. As the temperature of the hot wall is much higher in C3R1 case than in C1R1 case, the radiative transfer is much stronger in this case (see q R * w in Tab. 4.4) and the opposite effects of gas-gas and gas-wall radiative interactions shown in this latter case are amplified. The same analysis as for the case C1R1 can be achieved (see Sec. 4.4.1) and the results are similar, with an amplification effect for the conductive fluxes, as given in Tab. 4.4. Consequently, the averaged temperature is much higher than in C3 case. Figure 4.12 b shows that the usual log-law is not valid when radiation is accounted for. Moreover, no temperature log-law even appears between y + = 30 and y + = 200.

Influence of the Reynolds number

The results of the C4 and C4R1 cases, characterized by a higher Reynolds number than previously (see Tab. 4.1), are here compared with cases C1 and C1R1 in order to study the influence of the Reynolds number. The temperature profiles are plotted in Fig. 4.13. The curve T + (y + ) for case C4R1 in Fig. 4.13 b lies between the ones where radiation is not accounted for and the one of case C1R1. Increasing the Reynolds number between cases C1R1 and C4R1 has then moved the obtained wall-law closer to the usual one. There are two different explanations for this behavior. First, by increasing the Reynolds number, turbulent transport has been enhanced and its weight compared to the other energy transfer mechanisms is increased. This is shown in Fig. 4.14 a where turbulent transport is twice larger in case C4R1 than in case C1R1. One effect of the Reynolds number is therefore to relax the curve T + (y + ) obtained when radiation effects are considered towards the usual law of the wall which should be retrieved when radiative energy transfer is negligible.

In addition to this first effect, there is a second one that deals with a modification of the radiative energy transfer. In this specific configuration, the gas-gas and gas-wall radiative contributions to the total radiative power field in case C4R1 are modified as shown in in case C1R1. Magnitude of the gas-gas contribution decreases for the larger Reynolds number case while the magnitude of the gas-wall contribution increases. Since a larger Reynolds number leads to a larger gradient of the mean temperature close to the wall, the gas temperature at a given position gets closer to the average one in the core of the channel when the Reynolds number increases, while its difference with the near wall temperature gets larger. This explains the observed trends for the gas-gas and gas-wall radiative contributions. Consequently, the gas-gas effects on the wall conductive heat flux and the wall law dwindle while the gas-wall effects grow. This is seen in Tab. 4.4 where the conductive heat flux increase on the cold wall due to radiation is less in case C4R1 than in case C1R1. On the hot wall, the conductive heat flux even decreases between cases C4 and C4R1, showing that gas-wall effects overcome gas-gas effects there. Finally, the modification of the gas-gas and gas-wall contributions make the observed wall law go up compared to case C1R1. Among these two effects of the Reynolds number, the latter one is dominant in the present configuration where the radiative flux is much larger than turbulent transport (see Fig. 4.14). Would radiation remain dominant, increasing the Reynolds number even more could move the curve T + (y + ) above the usual wall law. This could not be check here due to the limitations of DNS on computational ressources with increasing Reynolds numbers.

Conclusion

In practical conditions of coupling between turbulent convection and radiation, there is no simple way for accurately predicting the averaged temperature profile and wall conductive flux, without undertaking a complete coupled computation. When the radiative flux within a turbulent gaseous medium is of the same order of magnitude or higher than the averaged turbulent convective flux, or the wall conductive flux, strong coupling effects occur within all the medium and at the walls. It was shown that the usual temperature profile and its corresponding log-law are generally no more valid within a turbulent boundary layer. The couplings between conduction, turbulent convection and radiation are complex and first strongly depend on both gas-gas and gas-wall radiative interactions. The global radiation effects are not easily predictable as gas-gas and gas-wall radiative interactions bring contributions of opposite sign to the wall conductive flux. Moreover, depending on the transversal optical thickness of the gaseous medium, complex effects, that are linked to multiple wall reflections, can also occur and strongly modify both the averaged temperature profiles and the conductive wall fluxes. Finally, the averaged temperature profiles and wall conductive fluxes were shown to also strongly depend on: i) the wall emissivities, that rule the wall reflection effects; ii) the temperature level, that controls the non linear radiative fluxes; iii) the Reynolds number that controls the weight of turbulence transport and the balance of gas-gas and gas-wall radiative effects. The development of simpler models for the determination of the temperature profile in the turbulent boundary layer is necessary.

Acknowledgements

The first author is supported by a doctoral fellowship from China Scholarship Council (CSC). This work was granted access to the HPC resources of [CCRT/CINES/IDRIS] under the allocation 2012-2b0164 made by GENCI (Grand Equipement National de Calcul Intensif). We also acknowledge Dr.

Vincent Moureau and the SUCCESS scientific group for providing the code YALES2 and support. Analysis of direct numerical simulations coupled to radiation in turbulent channel flow

Abstract

In the present paper, the role of radiative energy transfer on turbulent boundary layers is carefully analyzed, focusing on the effect on temperature fluctuations and heat turbulent transport. The work is based on direct numerical simulations of channel flows with hot and cold walls coupled to a Monte-Carlo method to compute the field of radiation power. In the studied conditions, the structure of the boundary layers are strongly modified by radiation. Temperature fluctuations and the turbulent heat flux are reduced, and new radiative terms appear in their respective balance equations. It is shown that they always counteract turbulence production terms. These effects are analyzed under different conditions of Reynolds numbers and wall temperatures where it is demonstrated that collapsing of wall-scaled profiles breaks down when radiation is considered. This is corrected by the introduction of a radiation-based scaling. Finally, the significant impact of radiation on heat turbulent transport is studied in terms of turbulent Prandtl number and a model for this quantity is developed based on the new proposed scaling and validated.

Introduction

Radiation plays an important role in many industrial applications, particularly combustion systems such as boilers, gas turbines, rocket engines and furnaces. For instance, in gas turbines, a crucial portion of the heat transferred from hot gas to the combustor solid walls comes from radiative energy transfer [START_REF] Lefebvre | Gas Turbine Combustion[END_REF]. The importance of radiation is even higher in modern gas turbines as the pressure ratio increases, which poses a severe problem on the cooling of the combustor walls. Moreover, radiation can influence the temperature distribution and hence the emission of pollutant in combustion systems. Therefore, an accurate prediction of radiation effect is important for the design of combustors. Among the studies of radiation effects in turbulent flows, great attention has been given to the interaction between turbulence and radiation (TRI). Two aspects of TRI can be identified: the effect of radiation on the temperature and species concentrations and vice versa. A comprehensive review about TRI is available in [START_REF] Coelho | Numerical simulation of the interaction between turbulence and radiation in reactive flows[END_REF][START_REF] Coelho | Turbulence-Radiation Interaction: From Theory to Application in Numerical Simulations[END_REF]. Regarding the effects of turbulence on radiation, it is observed that turbulence leads to an increase in the medium transmissivity [START_REF] Jeng | Radiative heat fluxes near turbulent buoyant methane diffusion flames[END_REF][START_REF] Gore | Spectral and total radiation properties of turbulent carbon monoxide/air diffusion flames[END_REF], the radiative power [START_REF] Coelho | Detailed numerical simulation of radiative transfer in a nonluminous turbulent jet diffusion flame[END_REF][START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF]) and heat loss [START_REF] Li | Importance of turbulence radiation interactions in turbulent diffusion jet flames[END_REF][START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF]). Coelho et al. (2003) reported that, in a non-luminous turbulent jet diffusion flame, TRI enhanced the heat losses "thesis_yufang" -2013/11/27 -19:31 -page 111 -#111

✐ ✐ ✐ ✐ ✐ ✐
Chapter 5 -DNS of turbulent channel flow: effect of radiation on second order moments 111 by a factor of 30% while a similar change in a luminous turbulent flame was also revealed in Ref. [START_REF] Tessé | Monte carlo modeling of radiative transfer in a turbulent sooty flame[END_REF]). Moreover, individual contributions to emission and absorption TRI have been isolated and quantified in a 1D premixed combustion system [START_REF] Wu | Direct numerical simulation of turbulence/radiation interaction in premixed combustion systems[END_REF], a homogeneous isotropic nonpremixed combustion system [START_REF] Deshmukh | Direct numerical simulation of turbulence-radiation interactions in homogeneous nonpremixed combustion systems[END_REF]) and a 1D turbulent nonpremixed flame [START_REF] Deshmukh | Direct numerical simulation of turbulence-radiation interactions in a statistically onedimensional nonpremixed system[END_REF].

By contrast to the former, only a few studies have been devoted to the effect of radiation on turbulence. Among them, Soufiani [START_REF] Soufiani | Temperature turbulence spectrum for high-temperature radiating gases[END_REF]) carried out a theoretical analysis of the influence of radiation on thermal turbulence spectra and it was concluded that radiation acted as a dissipation term and it could smooth the intensity of temperature fluctuations and modify the structure of the temperature variance spectrum. [START_REF] Damien | Analysis of the interaction between turbulent combustion and thermal radiation using unsteady coupled LES/DOM simulations[END_REF] also reported that radiation modifies the level of temperature fluctuations and homogenizes the spectral distribution of energy. Moreover, it was reported that the Reynolds stress and turbulence structure in supersonic shear layers were modified by radiation [START_REF] Ghosh | Effects of radiative heat transfer on the structure of turbulent supersonic channel flow[END_REF]. The objective of this paper is to investigate the effects of radiation in the boundary layer structure of turbulent channel flows. It has already been reported in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]) that radiation can significantly modify the mean temperature profile and consequently, the temperature wall law and the wall conductive heat flux. The different observed effects on the mean temperature profile have been understood thanks to the decomposition of radiation into wall-gas and gas-gas contribution. In this paper, Direct Numerical Simulations (DNS) of channel flows coupled with a reciprocal Monte Carlo method to deal with radiation from Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]) are analyzed focussing on radiation effects on higher-order statistical moments such as turbulent transport heat flux, enthalpy root-mean-square (RMS) and of their budget equations. After a detailed description of the studied problem in section 5.3, effects of radiation are analyzed in a first channel flow configuration in section 5.4.1. Then, changes of radiation effects with wall temperature difference and bulk Reynolds number are reported in section 5.4.2 where a new turbulent scaling is proposed. Finally, the associated results for the turbulent Prandtl number are shown in section 5.4.3 and a model based on the proposed scaling is derived.

Problem description

In order to study accurately the effects of radiation on the structure of turbulent boundary layers, direct numerical simulations of a planar channel flow coupled with a reciprocity Monte-Carlo method for radiation calculations have been considered. The set of governing equations in the fluid writes

∂ρ ∂t + ∂(ρu i ) ∂x i = 0, (5.1)
"thesis_yufang" -2013/11/27 -19:31 -page 112 -#112

✐ ✐ ✐ ✐ ✐ ✐ 112 5.3. PROBLEM DESCRIPTION ∂(ρu i ) ∂t + ∂(ρu i u j ) ∂x j = - ∂p ∂x i + ∂τ ij ∂x j + S i , (5.2) ∂(ρh) ∂t + ∂(ρu j h) ∂x j = ∂p ∂t - ∂q cd j ∂x j + P R (5.3) p = ρ r T, (5.4)
where ρ, u i , h, p and T are the fluid mass density, velocity components, enthalpy, pressure and temperature, respectively. h is expressed from the mixture thermal capacity at constant pressure c p : h = ∆h 0 + T T 0 c p (T ′ ) dT ′ , where T 0 is a reference temperature and ∆h 0 the corresponding standard formation enthalpy. The viscous shear stress tensor τ ij and the conductive flux vector q cd i write

τ ij = µ ∂u i ∂x j + ∂u j ∂u i - 2µ 3 ∂u k ∂x k δ ij , (5.5)
q cd i = -λ ∂T ∂x i , (5.6)
where µ is the dynamic viscosity, function of temperature computed like the mixture thermal capacity c p by the CHEMKIN package [START_REF] Kee | A fortran computer code package for the evaluation of gas-phase, multicomponent transport properties[END_REF][START_REF] Kee | CHEMKIN-II: A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase Chemical Kinetics[END_REF]. λ is the thermal conductivity, and is computed from the Prandtl number Pr=0.71. P R is the radiative power per unit volume. S i is a uniform forcing source term which acts as a pressure gradient term and drives the channel flow to obtain the desired bulk Reynolds number Re b . The set of governing equations is solved with the finite-volume solver YALES2 (Moureau et al. 2011b;Moureau et al. 2011c) under a low Mach-number approximation. As detailed in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF], the numerical setup is composed of a centered fourth-order spatial discretization and a fourth-order time integration. The computation of the radiation power is handled by an Optimized Emission-based Reciprocity Monte-carlo method (OERM) [START_REF] Zhang | Optimized emission-based reciprocity monte carlo method to speed up computation in complex systems[END_REF]. Gas radiative properties are calculated by using the weak absorption limit of CK model [START_REF] Soufiani | High temperature gas radiative property parameters of statistical narrow-band model for H 2 O, CO 2 and CO, and correlated-K model for H 2 O and CO 2[END_REF] in the studied cases that are characterized by high pressure.

The studied configuration, a fully developed turbulent channel flow with two isothermal walls, is shown in Fig. 5.1. The medium is a non-reacting CO 2 -H 2 O-N 2 gas mixture with a corresponding molar fraction of 0.116-0.155-0.729. Three computational cases from Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF], called here A, B and C, are defined in Tab. 5.1 by a set of bulk Reynolds number, pressure and wall temperatures (T w,c and T w,h ). When radiation is considered, theses cases are referred as A_RAD, B_RAD and C_RAD, respectively and the emissivity ε of the opaque walls is set to 0.8.

Results

Results for reference cases A and A_RAD

Results are first presented for cases A and A_RAD. Effects of radiation on the boundary layer structure are analyzed in terms of effects on the mean temperature field, on the enthalpy fluctuations and on the turbulent transport heat flux. Then, influence of temperature fluctuations on the radiative power is studied.

Mean temperature field

Owing to the small variation of mass density, mean velocity profiles (not shown here) are not affected by radiation when comparing cases A and A_R that are both characterized by a small wall temperature difference. However, the mean temperature profile, shown in Fig. 5.2 (a), is significantly modified by radiation over the whole domain. The mean temperature in wall units T + is defined as,

T + = |T -T w | T τ with T τ = |q cd w | ρ w c pw u τ , u τ = τ w ρ w 1/2 , (5.7)
where q cd w , ρ w , c pw and τ w are the mean conductive heat flux, mass density, thermal capacity and shear stress at the wall, respectively. The T + profiles of cases A and A_RAD are plotted in Fig. 5.2 (b) as functions of the normalized wall distance y + that writes,

y + = ρ w y u τ µ w , (5.8)
where µ w is the mean dynamic viscosity at the wall. The obtained temperature wall law for case A_R strongly deviates from that of case A, showing a significant effect of radiation on the thermal boundary layer structure in the considered conditions. A detailed analysis of such radiation effects on the mean temperature field in different channel flow conditions is available in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]. Opposite effects of radiation on the wall conductive heat flux and on the temperature wall law have been observed and they have been understood thanks to a decomposition of radiation into gas-gas and gas-wall contributions.

Fluctuations of enthalpy

Profiles of the enthalpy Favre root-mean-square, h rms = h ′′ h ′′ , in cases A and A_RAD are shown in Fig. 5.3 (a), where h rms is scaled by the center temperature defined as T c = (T w,c + T w,h )/2 and the center thermal capacity c p,c , also defined as an average between wall thermal capacities. In case A, without radiation, peaks of variance are located in the near wall regions as expected from standard boundary layer theory where production of turbulent fluctuations is maximal within the buffer layer. Because of the specific configuration where wall temperatures are different, a larger peak in h rms appears in the core of the channel where, as explained in Ref. [START_REF] Debusschere | Turbulent scalar transport mechanisms in plane channel and couette flows[END_REF], fluid pockets of high and low temperature converge from the hot and cold walls, respectively. In case A_R where radiation is accounted for, a significant reduction in enthalpy fluctuations is observed in the near wall region, especially on the hot side, and the central peak vanishes. Wall-scaled profiles of enthalpy rms against y + are presented in Fig. 5.3 (b), where the non-dimensional enthalpy rms h + rms is defined as

h + rms = h rms c pw T τ
(5.9)

Despite changes in the wall conductive heat flux are included in the wall scaling formulation through the definition of T τ , a large difference between results of cases A and A_R remains, indicating that the strong effect of radiation on fluctuations of enthalpy and temperature is a real modification of the boundary layer structure. This point is further investigated by analyzing the balance of the enthalpy variance transport equation which, in the studied configuration, writes

- ∂ ∂y (q cd ′ y h ′ ) - 1 2 ∂ ∂y (ρ v ′′ h ′′ h ′′ ) -ρ v ′′ h ′′ d h dy + q cd ′ i ∂h ′ ∂x i + h ′′ P R ′ + h ′′ ∂ ∂y (ρ v ′′ h ′′ ) = 0
(5.10)

where the terms on the left hand side are molecular diffusion, turbulent diffusion, production, molecular dissipation, correlation between enthalpy and radiative power fluctuations and a term proportional to h ′′ related to enthalpydensity correlation. These different terms are scaled by q cd w 2 /µ w and compared in Fig. 5.4 (a) and (b) for cases A and A_R (only the results on the cold side are shown since the hot side is similar). On the one hand, in case A, production and molecular dissipation terms are dominant and decrease away from the wall as expected in such standard conditions. On the other hand, in case A_RAD, a third dominant term appears in the balance of enthalpy variance in addition to the latter two, that is the enthalpy-radiative power correlation. Since this term appears as a negative contribution to the budget, it will be referred as radiative dissipation in the following. Hence, equilibrium between production and molecular dissipation away from the buffer layer for y + > 30 is replaced by a balance of production with molecular and radiative dissipations in the case with radiation. In the studied case, this equilibrium takes place sooner for y + > 20 and molecular dissipation remains weaker than radiative dissipation for y + > 50 approximatively. For both cases A and A_RAD, the term related to h ′′ , mean of mass-weighted fluctuating enthalpy, is negligible because of the small density variations. Regarding the scaled production, modification by radiation of the mean temperature field and of the turbulent heat flux (shown later) results in a reduction of this term in case A_R. This decrease in production and the presence of an additional radiation-related dissipative term in the budget of enthalpy variance explain the smaller level of enthalpy fluctuations shown in Fig. 5.3 for the case with radiation.

Wall-normal turbulent heat flux

Profiles of wall-normal turbulent heat flux ρ v ′′ h ′′ in cases A and A_R are presented in Fig. 5.5, where the wall-scaled turbulent heat flux is defined as

ρv ′′ h ′′ + = ρ v ′′ h ′′ ρ w u τ c pw T τ = ρ v ′′ h ′′ |q cd w |
.

(5.11)

The wall-normal turbulent heat flux is seen to decrease when radiation is accounted for. The aforementioned reduction in enthalpy fluctuations due to radiation is a first explanation why the turbulent transfer is less efficient. This effect can be filtered out by looking at the correlation coefficient between enthalpy and wall-normal velocity fluctuations shown in Fig. 5.6. It is observed that accounting for the local rms of enthalpy corrects the asymmetry of case A_RAD in Fig. 5.5 (a) and that correlation between enthalpy and wall-normal velocity is stronger with radiation. However, the correlation remains different in both cases, indicating another source of disagreement between the two cases than just the change in enthalpy rms.

The second reason is the requirement for the turbulent heat flux to fulfill the mean balance equation for energy which writes,

- ∂ ∂y q cd y + ρ v ′′ h ′′ + P R = 0.
(5.12)

Introducing the radiative flux vector q R i , the mean energy balance equation writes as a constant sum of energy fluxes, where q cd w and q R w are the wall conductive heat flux and wall radiative flux, respectively. Hence, the turbulent heat flux is constrained to equilibrate with conductive and radiative fluxes, establishing a strong two-way coupling between the mean and the fluctuations of the temperature and radiative power fields. More information about the profiles of energy fluxes is available in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF].

q cd y + ρ v ′′ h ′′ + q R y = q cd w + q R w , (
The balance equation of the turbulent heat flux ρ v ′′ h ′′ writes

∂ ∂y τ ′ 22 h ′ -q cd ′ y v ′ - ∂ ∂y ρ v ′′ v ′′ h ′′ -ρ v ′′ v ′′ d h dy -τ ′ 2i ∂h ′ ∂x i -q cd ′ i ∂v ′ ∂x i +v ′′ P R ′ + p ′ ∂h ′ ∂y - ∂p ′ h ′ ∂y + h ′′ ∂ ∂y ρ v ′′ v ′′ + v ′′ ∂ ∂y ρ v ′′ h ′′ = 0.
(5.14)

The seven terms on the left hand side are molecular diffusion, turbulent diffusion, production, molecular dissipation, correlation between wall-normal velocity and radiative power fluctuations, pressure contribution and a term related to the average of Favre fluctuations of enthalpy and velocity. These terms, except for the last one which is again negligible, are shown for cases A and A_RAD in Fig. 5.7 (a) and (b) where they are scaled by |q cd w | τ w /µ w . In case A, predominant terms are production and the pressure contribution and molecular dissipation accounts for the remaining balance. As shown in Ref. [START_REF] Kasagi | Direct numerical simulation of passive scalar field in a turbulent channel flow[END_REF]), the two terms in the pressure contribution are of the same order of magnitude. When radiation is taken into account, most of the terms are reduced in magnitude and broader. For y + > 25, the balance is split into production, the pressure contribution, the wall-normal velocity-radiative power correlation term and molecular dissipation which has the smallest contribution among these four terms. The pressure term is dominant for y + < 100 until the velocity-radiative power correlation term takes over.

Radiative power field

The mean profile of the radiative power is shown in Fig. 5.8 (a). It is positive (resp. negative) in the very near wall region on the cold side (resp. hot side). Further away from the wall, the radiative power changes sign twice around y + = 20 and 100 on both sides. This shape of the mean radiative power has been understood thanks to a separation of gas-gas and gas-wall contributions to radiation done in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]. One of the most studied aspect of turbulence-radiation interaction is the effect of turbulent fluctuations on the mean radiative power [START_REF] Coelho | Turbulence-Radiation Interaction: From Theory to Application in Numerical Simulations[END_REF]. Indeed, the mean radiative power writes where the spectral intensity I ν , the equilibrium spectral intensity I • ν and the spectral absorption coefficient κ ν are integrated over frequency ν and solid angle Ω ranges. Equation 5.15 can alternatively be expressed as

P R = +∞ 0 4π κ ν I ν dΩ -4πκ ν I • ν (T ) dν , ( 5
P R = +∞ 0 4π κ ν I ν + κ ′ ν I ′ ν dΩ -4π κ ν I • ν (T ) + κ ′ ν I • ′ ν (T ) dν , (5.16)
Hence, local and non-local correlations between spectral intensities and absorption coefficients appear in the expression of the mean radiative power. Besides, radiation is highly non-linear with temperature. For all theses reasons, in most cases, P R cannot be computed from the mean temperature field, i.e.

P R ({T }) = P R ({T }).

(5.17)

To investigate this effect, the radiative power computed from the mean temperature field is also shown in Fig. 5.8 (a). No difference is seen with the exact profile, showing that turbulent fluctuations of temperature are not intense enough to perturb the mean radiative power. This result can be generalized to most thermal boundary layers when the wall and bulk temperature ratio is low. Indeed, it is demonstrated in appendix that, in air, the level of temperature fluctuations is very small when with a small temperature ratio and it even decreases with Reynolds number. However, when the temperature ratio increases, the fluctuation level can be as high as 30%, which might be high enough to have an effect on the mean radiative power. Moreover, such an interaction between turbulence and radiation is mostly observed in combustion applications [START_REF] Coelho | Turbulence-Radiation Interaction: From Theory to Application in Numerical Simulations[END_REF] where the heat release stemmed from chemical reactions greatly enhances temperature fluctuations. Although the mean radiative power is not influenced by temperature fluctuations in the present case, fluctuations of the radiative power field modify the balance of enthalpy variance and turbulent heat flux transport equations, Eqs. 5.10 and 5.14, and hence interact with the mean temperature field. In order to compute the relevant radiative power root-mean-square P R rms , it is necessary to subtract the standard error of the Monte-Carlo method from the total rms data. Indeed, since the Monte-Carlo method is a statistical approach, the computed time-averaged variance of radiative power writes (5.18) where the instantaneous radiative power P R = P R phys + P R err is composed of the real physical value of radiative power, P R phys , that is estimated by the Monte-Carlo method with a controlled error, and of a statistical error P R err . Assuming independency, between physical and statistical fluctuations, the root-meansquare of radiative power is computed as using the statistical estimation of P R 2 err provided by the Monte-Carlo approach. The profile of radiative power root-mean-square P R rms is plotted in Fig. 5.8 (b). Contrary to the mean radiative power that abruptly vanishes away from the walls, fluctuations of radiative power do not disappear and their magnitude presents only slight variations for y+ > 20.

P R ′2 = P R ′2 phys + P R 2 err + 2P R ′ phys P R err ,
(P R rms ) 2 = P R ′2 -P R 2 err , (5.19)
As shown previously, the fluctuating radiative power introduces a new radiative dissipation term in Eq. 5.10 that writes h ′′ P R ′ . Since radiation power is the difference between absorbed and emitted powers, P R = P R a -P R e , the radiative dissipation term is also composed of two terms:

h ′′ P R ′ = h ′′ P R ′ a -h ′′ P R ′ e .
(5.20)

As shown in Ref. [START_REF] Ammouri | Effects of temperature and concentrations in turbulent gas flows on combined radiative and conductive wall fluxes[END_REF], the emitted part always leads to a negative contribution to the total radiative dissipation term because a positive (resp. negative) fluctuation of enthalpy yields a higher (resp. smaller) emitted radiative power, while the absorbed part is mainly positive. The socalled radiative dissipation could then be positive in some conditions. However, the corresponding correlation coefficient between enthalpy and radiative power shown in Fig. 5.9 (a) is negative and indicates that fluctuations of radiative power are fairly correlated with local fluctuations of enthalpy or temperature, especially close to the wall. The radiative dissipation is then here dominated by its emission part.

In the balance equation of the turbulent heat flux (Eq. 5.14), the term v ′′ P R ′ related to the correlation coefficient shown in Fig. 5.9 (b) appears if radiation is taken into account. In the present case, this term is positive in the whole computational domain although there is no reason to generalize this result. It is worth noting that, by introducing the radiative energy flux, this term can also be split into two terms as done for the conductive heat flux in Eq. 5.14:

v ′′ P R ′ = - ∂ ∂x i (q R ′ i v ′ ) + q R ′ i ∂v ′ ∂x i (5.21)
It seems interesting to understand the sign of the term v ′′ P R ′ in view of the fair correlation between enthalpy and radiative power fluctuations. Indeed, one could argue that the fluctuation of wall-normal velocity leads to a fluctuation of enthalpy due to convection which finally generate a variation in the radiative power. In order to verify this effect, a first attempt could consist in passing over the convected enthalpy fluctuations and directly approximate the variation of radiative power as a result from convection and the shape of the mean profile P R . Turbulent diffusion of the mean radiative power with a gradient assumption would then lead to

v ′′ P R ′ ∼ - ∂P R ∂x i (5.22)
Although this approximation gives the right sign for the term v ′′ P R ′ in the very near wall region on both sides, it becomes rapidly erroneous for y + > 25 where the monotonicity of the mean radiative power profile changes. Skipping the enthalpy fluctuations to understand the correlation between wall-normal velocity and the radiative power is therefore not appropriate. Another approach is proposed by considering the fourth-order cross-moment v ′′ h ′′ h ′′ P R ′ which introduces enthalpy fluctuations as the vessel that could explain the correlation mechanism between velocity and radiative power variations. Two high-order correlation coefficient C 1 and C 2 are then defined as

v ′′ h ′′ h ′′ P R ′ = C 1 v ′′ h ′′ h ′′ P R ′ (5.23) v ′′ h ′′ h ′′ P R ′ = C 2 h ′′ h ′′ v ′′ P R ′ (5.24)
Finally, the term v ′′ P R ′ can be indirectly related to the high-order statistical moment by writing

v ′′ P R ′ = C v ′′ h ′′ h ′′ P R ′ h ′′ h ′′ , (5.25)
where C, defined as C 1 /C 2 , is computed from Eq. 5.25 and plotted in Fig. 5.10. The profile of the coefficient C is regular, positive and varies from approximatively unity close to the walls to 2.7 in the core of the channel. Since the radiative dissipation term h ′′ P R ′ is negative, it can then be postulated from Eq. 5.25 that the term v ′′ P R ′ appearing in Eq. 5.14 may always have the opposite sign of the turbulent flux ρ v ′′ h ′′ . In the present configuration, this explains the positive sign of the correlation between wall-normal velocity and radiative power.

Being of opposite sign with the turbulent heat flux, the term v ′′ P R ′ would therefore always be a loss term in Eq 5.14. This loss mechanism described by Eq. 5.25 is as follows: A positive or negative variation of wall-normal velocity generates a fluctuation of enthalpy whose sign is determined by the direction of the turbulent heat flux; And this enthalpy fluctuation then induces a variation in radiative power of opposite sign.

Summary of radiation effects in turbulent channel flow

The coupling between the temperature and radiative fields has been carefully detailed through the analysis of DNS/Monte-Carlo results obtained in cases A and A_R. Figure 5.11 sums up the different interactions between mean (Favre or Reynolds average) and fluctuating fields of velocity, temperature and radiation observed in the present channel flow simulations. The feed back of mean temperature on velocity is done through the variation of temperature-dependent properties such as mean density and dynamic viscosity. In cases A and A_R characterized by a relative small difference of wall temperatures, this effect is negligible although it is present in cases with larger temperature difference as reported in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]). On the other hand, impact of temperature fluctuations on the velocity field mainly through density variations is negligible in all investigated cases, making Reynolds and Favre averaging operations similar. In the specific configuration of a channel flow, convection by the mean velocity field does not influence the temperature field directly. It has obviously a major role in any general flow configuration. However, at least through the turbulent heat flux and the production of enthalpy variance, the fluctuations of velocity play a critical role in the profile of mean and root-mean-square of temperature. As outlined in the mean radiative power strongly depend on each other. The aforementioned modification of production in Eq. 5.10 is an indirect effect of radiation involving the two-way interaction between P R , T and T ′ . Another indirect effect of radiation on the mean temperature comes from the interaction between fluctuations of temperature and radiation. Accounting for radiation makes additional loss terms appear in the balance equations of enthalpy variance and turbulent heat flux, which makes turbulent transport of the mean temperature field less efficient. The variations of radiative fields are also determined from the level of mean temperature and mean radiative fields. Finally, it has been observed that the influence of fluctuations of the temperature and radiative fields is not noticeable on the mean radiative field.

Effects of wall temperature difference and Reynolds number

In this section, cases A, B and C are compared with and without accounting for radiation to understand the effects of wall temperature difference and Reynolds number. Analysis of the results focuses on the turbulent heat flux, the enthalpy root-mean-square and the main contributing terms in their balance equations.

A new scaling is proposed to improve collapsing of several profiles.

Effect of wall temperature difference

In comparison to the cases A and A_R, cases B and B_R are characterized by a large wall temperature difference (see Tab. 5.1). Mean temperature profiles of these four cases are shown in Fig. 5.12. Small differences between cases without radiation are due to noticeable effects of mean density variations. When ac- counting for radiation, the modification of the mean temperature in both T and T + profiles is stronger between case B and B_R because of the amplification of radiation at high temperature.

Profiles of mean velocity reported in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]) exhibit an more significant impact of the mean density on cases B and B_R. The turbulent kinetic energy in wall units k + = 0.5 u ′′ i u ′′ i /u 2 τ is plotted in Fig. 5.13 (a). Similarly, the variations of mean density prevent collapsing of the profiles, especially when comparing cases with small and large wall temperature difference. The effect of radiation through the change in mean temperature profile and, hence, mean density is small. In order to improve the agreement between profiles of different cases, an alternative semi-local scaling is applied where non-uniformity of gas properties is considered by using their local values (Huang et al. 1995;[START_REF] Coleman | A numerical study of turbulent supersonic isothermal-wall channel flow[END_REF][START_REF] Dailey | Large eddy simulation of constant heat flux turbulent channel flow with property variations: Quasideveloped model and mean flow results[END_REF]. The corresponding friction velocity u * τ and friction temperature T * τ are then defined as: filtered out and to validate whether the same physical interpretation can be use to understand the shape of the profiles. Here, it confirms that radiation affects turbulent kinetic energy indirectly through the change of temperaturedependent gas properties. Henceforward, semi-local scaling is used to normalize turbulent quantities in replacement of standard wall-scaling. Using semi-local scaling, the scaled turbulent heat flux, perature and non-linear, which explains why the deviation from cases without radiation is larger in case B_R because of its higher temperatures.

u * τ = τ w ρ 1/2 , T * τ = |q cd w | ρ c p u * τ , ( 5
ρv ′′ h ′′ * = ρ v ′′ h ′′ ρu * τ c p T * τ = ρ v ′′ h ′′
The three dominant terms in the balance equation (Eq. 5.14) of the turbulent heat flux, i.e. production, pressure and radiation terms, are also reported in Fig. 5.14. Following the semi-local scaling approach, these quantities are scaled by |q cd w |τ w /µ. The collapsing of production and pressure-term profiles corresponding to cases A and B is good except for the hot side in case B which presents a small difference that could not be corrected entirely by semi-local scaling.

In cases with radiation, the scaled production that is related to

-v ′′ v ′′ * d h * dy *
is smaller in magnitude. Since, as shown for the turbulent kinetic energy k, the semi-local scaled root-mean-square of wall-normal velicity v ′′ v ′′ * is barely affected by radiation through the mean density variation, the modification of the production magnitude is due to the change of scaled mean enthalpy gradient d h * dy * which is tightly connected to dT + dy + . Looking at the gradient of the scaled mean temperature T + in cases with radiation for y + between 10 and 30 in Fig. 5.12 (b), where the production peaks, the same variations between different cases are almost retrieved for the peak value of scaled production in Fig. 5.14 (b). This demonstrates the major role of the wall-scaled enthalpy gradient in this term.

As for the pressure term in Fig. 5.14 (c), since it mostly compensates production for y * < 100, the same trend is retrieved. The scaled correlation between wall-normal velocity and radiative power that appears in Eq. 5.14 is shown in Fig. 5.14 (d). For both A_R and B_R, this term does not vary much beyond y * = 50 and is larger in case A_R. The explanation is not obvious as the magnitude of this radiative loss term is determined by multiple phenomena as described in section 5.4.1.4.

Radiation-based scaling

It has been shown that semi-local scaling is working well for cases without radiation by accounting for the effects of variable flow properties. However, it is not able to handle that, under different conditions, radiation changes in magnitude and in nature (gas-gas versus gas-wall contributions [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]). Consequently, each case appears as different from the others although the same physical interpretation might stand to explain the shape and magnitude of the observed profiles as seen for the production term.

In order to derive a new scaling that could account for the variability of radiation effects, let's first consider the mean energy balance equation without radiation:

q cd y (Y ) + ρ v ′′ h ′′ (Y ) = q cd w .
(5.29)

Consequently, outside of the viscous sublayer, where the conductive heat flux is negligible, the turbulent heat flux is the same as the wall conductive flux q cd w .

" It is then the right quantity for scaling of the turbulent heat flux and to define the friction temperature. Similarly, when radiation is taken into account, the mean energy balance equation in Eq. 5.13 is written as (5.30) where the change of radiative flux in respect to the wall radiative flux ∆q R (Y ) is given by

q cd y (Y ) + ρ v ′′ h ′′ (Y ) = q cd w -∆q R (Y ) ≡ q † (Y ) ,
∆q R (Y ) = q R (Y ) -q R w .
(5.31)

The definition ∆q R (Y ) automatically suppress any wall-wall radiative contribution that, regarding the fluid, is a passive energy transfer. For cases involving radiation, the sum of conductive and turbulent heat flux q † (Y ) is a function of wall distance whose shape is determined by the intensity and nature of radiative energy transfer. For cases without radiation, q † (Y ) turns back to the wall conductive flux. Normalizing by |q † (Y )| the turbulent heat flux, the corresponding profiles of ρv ′′ h ′′ † are shown in Fig. 5.15 (a). In comparison to Fig. 5.14 (a), the agreement between the different cases is significantly improved. Thanks to the radiation-based scaling, the scaled turbulent heat flux represents the relative strength of turbulent transport compared to conduction for all cases. That is why all profiles now present the same shape and reach a plateau. The remaining discrepancies concern the level of this plateau that is related to the remaining contribution of the conductive heat flux which, in the studied configurations with hot and cold walls, does not tend to zero in the core of the channel. The level of the plateau is then determined by the mean temperature gradient at the center of the channel and is different for each case. It is expected that, for symmetrical channel flows with the same wall temperatures or external boundary layers, this effect is not present, making the proposed scaling even more efficient.

The three main terms in the balance equation of turbulent heat flux are also shown in Fig. 5.15, where production, pressure and radiation terms are now scaled by |q † (Y )|τ w /µ. The agreement between profiles is also much better when the radiation-based scaling is considered. In order to understand why, the scaled production term is written as has improved the collapsing of production profiles, its main counter-balancing terms that are the pressure and radiation terms benefit from the same effect although it is not striking for the latter. The enthalpy root-mean-square and the main terms in the enthalpy variance balance equation (Eq. 5.10), production, molecular and radiative dissipations, are shown in Fig. 5.16 using semi-local scaling (the transport equation terms are scaled by q cd w 2 /µ). Once again, the good collapsing of profiles corresponding to cases without radiation demonstrate the efficiency of semi-local scaling to account for variable flow properties effects. Nonetheless, when radiation is taken into account, rms of enthalpy, its production and molecular dissipation are totally different. As shown previously, then deviation from cases without radiation is stronger for higher temperatures. The same results using radiation-based scaling are presented in Fig. 5.17 (5.33) the normalized enthalpy h † rms is scaled by c p T † τ while the balance equation terms are scaled by the quantity q † (Y ) 2 /µ. The radiation-based scaling allows similar shapes of h † rms to be retrieved for all cases. Cases with and without radiation seem to collapse together separately although case B_R deviates from case A_R for y * > 75. As the proposed radiation-based scaling mainly accounts for variations of the turbulent heat flux profiles, the fact that the enthalpy rms is not directly expressed as a function of the turbulent heat flux by a simple expression can explain this discrepancies. On the other hand, the scaled production that writes (5.34) is related to the turbulent heat flux and the obtained agreement between production profiles in Fig. 5.17 (b) is good. Production of enthalpy variance is mainly compensated by molecular and radiative dissipations which then individually show a fair agreement of their profiles using the new scaling. As radiation introduces a new significant loss term in the balance equation of enthalpy variance, the agreement between profiles of molecular dissipation seems only valid for cases without and with radiation separately. However, switching progressively from a case without radiation to cases with more and more significant radiative effects, one would expect a continuous transition. This trend is more visible in other cases and is detailed in the next section on Reynolds number effect. Finally, it should be mentioned that the use of the radiation-based friction temperature T † τ is not able to improve the agreement between mean temperature profiles (see Fig. 5.18). Like the enthalpy root-mean-square, the mean temperature cannot be directly expressed as a function of the turbulent heat flux. Introducing the turbulent thermal diffusivity, the turbulent heat flux is in fact related to the gradient of mean temperature. Theses quantities can alternatively be connected through the mean energy balance equation in Eq. 5.30. Therefore, although the radiation-based friction temperature may properly scale the gradient of mean temperature, its dependence on wall distance prevents it from appearing straightforwardly in the expression of the mean temperature. transport, which makes the temperature profile in case C_R closer to those of cases without radiation. On the other hand, a modification of the nature of radiative energy transfer has also been observed [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF].

-ρv ′′ v ′′ d h d y |q † (Y )|τ w /µ ≈ -ρ v ′′ v ′′ * Pr 1 -ρv ′′ h ′′ † , ( 5
T † τ = |q † (Y )| ρ c p u τ ,
-ρ v ′′ h ′′ d h d y q † (Y ) 2 /µ ≈ -ρv ′′ h ′′ † Pr 1 -ρv ′′ h ′′ † ,

Effect of Reynolds number

Regarding the profiles of mean turbulent heat flux and its balance equation main terms shown in Fig. 5.20 using semi-local scaling, good agreement of profiles without radiation is obtained while results of case C_R sit between profiles of cases without radiation and the ones of case B_R because of the weakened radiative effects at higher Reynolds number. The same profiles are reported in Fig. 5.21 using the radiation-based scaling. Agreement between profiles is improved although results for case C_R are still noticed to remain between case B_R and case A/B. Similar observations are obtained for the enthalpy root-mean-square and for its production and dissipation terms shown in Fig. 5.22 and 5.23 using semi-local scaling and radiation-based scaling, respectively. Contrary to the comparison of cases A_R and B_R, it clearly appears in Fig. 5.23 (d) for the scaled molecular dissipation that there is a transition between cases without radiation and cases with stronger and stronger radiative energy transfer. Consequently, although profiles of production correctly collapse with the proposed radiation-based scaling, its necessary splitting into two main terms, when radiation is taken into account, prevents collapsing for molecular and radiative dissipation profiles separately. The balance between molecular and radiative dissipations is controlled by the importance of radiation so that a smooth transition between profiles of different cases can be expected when radiative energy transfer increases from being negligible to strongly dominant. With the same arguments, the profiles of pressure and radiation terms in the balance equation of the turbulent heat flux are not supposed to entirely collapse and the balance between the two terms is determined by the weight of radiation effects.

Effects on the turbulent Prandtl number

Previously results have outlined the effect of radiation on the turbulent heat flux. In order to model turbulent heat transfer, the turbulent Prandtl number Pr t is a very important quantity that is extensively used to relate the thermal eddy diffusivity a t to the turbulent eddy viscosity ν t . It is defined as

Pr t = ν t a t = -ρ u ′′ v ′′ -ρ v ′′ T ′′ dT /dY du/dY . (5.35)
The definition could alternatively be based on the gradient of Favre average of temperature and velocity. As mentioned earlier, the difference between Reynolds and Favre averaging approaches is small in the studied conditions. The profiles of Pr t are presented in Fig. 5.24 for all studied cases. The profile of the turbulent Prandtl number in boundary layers has been widely studied [START_REF] Kays | Turbulent Prandtl number. Where are we?[END_REF]) and results for cases without radiation are consistent with channel flow simulations in the literature with air (Pr= 0.71) and fixed temperature at the wall (Kim and Moin 1987;[START_REF] Kong | Direct numerical simulation of turbulent thermal boundary layers[END_REF]: The value is close to unity at the wall and decreases along the wall-normal distance, reaching an asymptotic value for high Reynolds number flows. Deviation of case B from cases A and C is probably due to mean density effects. Nonetheless, cases with radiation exhibit a much stronger difference where the obtained profiles of turbulent Prandtl number are all above the standard ones without radiation. Profiles for cases A_R, B_R and C_R on different sides of the channel are all different from each other but have the same pattern: Pr t decreases from a value larger than unity at the wall until y * ≈ 15, then increases between y * = 15 and 75 to reach a plateau beyond. The value y * = 75 also corresponds to the value where the radiation-based turbulent heat flux reaches a plateau (see Figs. 5.14 and 5.21). This is not surprising since Pr t is related to ρv ′′ h ′′ † :

Pr t ∼ dT /dY -ρ v ′′ T ′′ ∼ 1 -ρv ′′ h ′′ † -ρv ′′ h ′′ †
(5.36)

Finally, the effects of Reynolds number and wall temperature difference are consistent with previous results: Case C_R is closer to results without radiation while the hot wall profile in case B_R presents the largest difference. In Reynolds Averaged Numerical Simulations (RANS), the turbulent Prandtl number is either taken as constant or calculated by a more comprehensive model. In the latter, the model usually uses the time scale ratio R = τ t /τ u , where the velocity turbulent time scale τ u and thermal turbulent time scale τ t are defined as

τ u = k ǫ , τ t = h ′′ h ′′ 2ǫ t , (5.37)
and ǫ and ǫ t are dissipation terms of the turbulent kinetic energy k and enthalpy variance, respectively. When radiation is accounted for, there are two different types of dissipation, molecular and radiative dissipation, so that three different thermal time scales could be defined: the first based on molecular dissipation only, the second based on radiative dissipation only and the third one based on the sum of theses two terms. The different corresponding time scale ratios are shown in Fig. 5.25 in case A and A_R. The common time scale ratio based on thermal molecular dissipation is larger in case A_R while the one based on radiative dissipation is much smaller. Finally, the combination of the two dissipative terms in the definition of R leads to the shortest thermal turbulent time scale. As radiation dissipation dominates, this latter definition is close to the one based on radiative dissipation only. Different models for the thermal turbulent diffusivity a t , which takes the time scale ratio R as an input and are then able to predict variable turbulent Prandtl number, are now considered. These models are completed with four transport equations to compute the different turbulent time scales. Here, DNS data are used to assess their accuracy. The turbulent diffusivity a t is modeled as

a t = C λ k τ m f λ , (5.38)
where τ m is a mixed time scale expressed from τ u and R, C λ a model constant usually set to retrieve the asymptotically constant turbulent Prandtl (∼ 0.85) in air at high Reynolds number and f λ is a damping function accounting for low-Reynolds effects very close to the wall. Two formulations for τ m are considered.

The first one (referred as model A) used in Refs. [START_REF] Nagano | A two-equation model for heat transport in wall turbulent shear flows[END_REF][START_REF] Soufiani | Radiation effects on turbulent heat transfer in channel flows of infrared active gases[END_REF] writes

τ m ∼ τ u R m with m = 0.5, (5.39)
while the second formulation (referred as model B) used in Ref. [START_REF] Abe | A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows-ii. thermal field calculations[END_REF]) can be written

τ m ∼ τ u R R + C m .
(5.40)

For each of these models, the thermal turbulent time scale can be evaluated either from the molecular dissipation only (modeled referred and A1 and B1) or from the total dissipation, sum of molecular and radiative dissipative terms (modeled referred as A2 and B2). The latter approach was already proposed in Ref. [START_REF] Ammouri | Effects of temperature and concentrations in turbulent gas flows on combined radiative and conductive wall fluxes[END_REF]. Models A1, A2, B1 and B2 are a priori tested using data from DNS to compute a t . Models constant values and damping "thesis_yufang" -2013/11/27 -19:31 -page 143 -#143

✐ ✐ ✐ ✐ ✐ ✐
Chapter 5 -DNS of turbulent channel flow: effect of radiation on second order moments 143 functions for A1 and A2 are taken from Ref. [START_REF] Soufiani | Radiation effects on turbulent heat transfer in channel flows of infrared active gases[END_REF]) while, for models B1 and B2, they are taken from Ref. [START_REF] Abe | A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows-ii. thermal field calculations[END_REF]. The predicted profiles of a t are combined with the turbulent viscosity profiles extracted from DNS results to compute profiles of turbulent Prandtl number Pr t that are given in Fig. 5.26.

In case A, without radiation, predictions of models A1 and B1 are similar and tend to an asymptotic value of Pr t . Agreement with the DNS data is correct although the turbulent Prandtl number obtained in DNS does not reach a plateau and keeps decreasing. This behavior can be attributed to low-Reynolds effects in the core of the channel but also to the specific configuration of different wall temperatures, which leads to an increase of enthalpy fluctuations, starting at y + ≈ 100, in the center of the channel (see Fig. 5.3).

With radiation, results of model B1 (and A1, not shown) shows that considering only the molecular dissipation significantly underestimates the turbulent Prandtl number and even predicts the wrong trend, the predicted values of Pr t being smaller. On the other hand, when both molecular and radiative dissipations are used in the definition of the thermal time scale, the turbulent Prandtl number predicted by models A2 and B2 is correctly boosted by radiation effects but the deviation from the DNS profile remains important.

Using previous results with the radiation-based scaling, another model for Pr t is proposed to account for radiation effects in turbulent boundary layers. In section 5.4.2, it has been demonstrated that the non-dimensional turbulent heat flux, when normalized with radiation-based scaling, agree well between different cases with and without radiation, i.e.,

ρv ′′ h ′′ † R ≈ ρv ′′ h ′′ † 0 , (5.41)
where subscript indices R and 0 correspond to quantities in cases with and without radiation, respectively. For cases with radiation, the thermal turbulent diffusivity can then be written

a t,R = 1 ρ R c p R ρv ′′ h ′′ R dT R /dy ≈ 1 ρ R c p R |q † R (Y )| dT R /dy ρv ′′ h ′′ † 0 (5.42)
For cases without radiation, ρv ′′ h ′′ † 0 is expressed as

ρv ′′ h ′′ † 0 = ρv ′′ h ′′ 0 |q cd w,0 | = a t,0 a 0 + a t,0 dT 0 /dy |dT 0 /dy| = ν t,0 a 0 Pr 0 t + ν t,0 dT 0 /dy |dT 0 /dy| , (5.43)
where a 0 is the thermal diffusivity. From Eqs. (5.42) and (5.43), the turbulent Prandtl number Pr t,R when radiation effects are taken into account can be calculated as Neglecting radiation effects on the velocity field, ν t,0 can then be replaced with ν t,R . Assuming also that a R ≈ a 0 , finally gives

Pr t,R = ν t,R a t,R = ρ R c p R |dT R /dy| |q † R (Y )| ν t,R a 0 Pr 0 t + ν t,0 ν t,0 . ( 5
Pr t,R = ρ R c p R |dT R /dy| |q † R (Y )| a R Pr 0 t + ν t,R , (5.45)
Effects of radiation are included through q † R (Y ) which is related to the radiative flux and through the change of the mean temperature profile in the term |dT R /dy|. When no radiative energy transfer is considered, the previous relation simplifies to Pr 0 t , profile of turbulent Prandtl number in cases without radiation. This last piece of the model can either be given by an algebraic formula to predicted by a more complex model. Here, the formula proposed in Ref. [START_REF] Kays | Turbulent Prandtl number. Where are we?[END_REF]) is considered:

Pr 0 t = 1 0.5882 + 0.228(ν t,R /ν R ) -0.0441(ν t,R /ν R ) 2 1 -exp -5.165 ν t,R /ν R . (5.46)
Equations 5.45 and 5.46 form a turbulent Prandtl model for turbulent boundary layers with radiation that can be used in RANS simulations but also in wallmodel for large eddy simulations. Using DNS data to evaluate Eq. 5.45, results of the proposed model are compared with DNS results in Fig. 5.27 for all studied cases. The agreement of the model with DNS is impressive, showing the importance of the radiation-based scaling.

Conclusion

Several direct numerical simulations of turbulent channel flows composed of burnt gases under different conditions (Reynolds number, wall temperature difference) with and without radiative energy transfer, computed with a Monte-Carlo method, are analyzed. Beside the previously observed modification of the mean temperature profiles and associated wall laws, fluctuations of temperature and the turbulent heat transport are significantly reduced by radiation effects in the studied conditions for two main reasons. On the one hand, radiation indirectly influence production terms that appear in the balance equation of these quantities by changing gradients of the mean temperature field. On the other hand, the additional radiative power source term in the energy balance equation creates new terms in transport equations of enthalpy variance and turbulent heat flux. The two new terms, correlations between radiative power and enthalpy fluctuations and between radiative power and wall-normal velocity fluctuations, appears as loss terms in these equations. There is therefore a strong interaction between the fluctuations of radiative power and of temperature, and the latter then influences the mean temperature field through turbulent transport. However, the widely studied impact of turbulent fluctuations on the mean radiative power is negligible in the considered channel flows where no combustion takes place. The impact of radiation on the turbulent heat flux, the enthalpy variance and their transport equation terms is shown to be different in each case, for different wall temperature and different Reynolds. The usual profiles collapsing of wallscaled variables under these different conditions then falls apart when radiation is taken into account. This effect is corrected by introducing a radiation-based scaling instead. Good agreement between profiles of scaled turbulent heat flux and production terms is obtained for all studied cases. The balance between scaled molecular dissipation (or pressure term) and radiation loss term, that both compensate production, remains determined by the importance of radiative energy transfer. Finally, using radiation-based scaling, a model for the turbulent Prandtl number in boundary layers with radiation effects is proposed and validated. This model can be used in RANS or wall-modeled LES when radiation is able to modify the boundary layer structure. In a turbulent thermal boundary layer, the peak of temperature rms is located in the buffer layer. The peak value and position have been shown to depends on both Reynolds number and Prandtl number (Kim and Moin 1987;[START_REF] Kawamura | DNS of turbulent heat transfer in channel flow with low to medium-high Prandtl number fluid[END_REF][START_REF] Kawamura | DNS of turbulent heat transfer in channel flow with respect to reynolds and prandtl number effects[END_REF]Abe et al. 2004). In air with Pr= 0.71, optimum value of temperature rms in wall units and its position become Reynolds-independent at sufficiently high Reynolds numbers [START_REF] Kasagi | Direct numerical simulation of passive scalar field in a turbulent channel flow[END_REF][START_REF] Kawamura | DNS of turbulent heat transfer in channel flow with low to medium-high Prandtl number fluid[END_REF][START_REF] Kawamura | DNS of turbulent heat transfer in channel flow with respect to reynolds and prandtl number effects[END_REF][START_REF] Kong | Direct numerical simulation of turbulent thermal boundary layers[END_REF]Abe et al. 2004). Therefore, the maximum temperature rms writes

T max rms = T rms (y + BL (Re, Pr)) = T +,max rms (Re, Pr) T τ , (5.47)
where y + BL is the position of the peak and T +,max rms is the maximum temperature rms, both in wall units. For air, it is found that y + BL ≈ 12 and T +,max rms ≈ 2.5. Common approximations of the friction coefficient and Nusselt number in ducts are

c f = τ w ρu 2 b ≈ αRe -0.2 , Nu = |q cd w |D h |T w -T b |λ ≈ βRe 0.8 Pr 0.5 (5.48)
where u b is the bulk velocity, T b the bulk temperature, D h the hydraulic diameter of the duct and Re= ρu b D h /µ. The constant coefficients are α = 0.023 and β = 0.022 in air. The friction temperature can then be expressed as

T τ = α -0.5 βRe -0.1 Pr -0.5 |T w -T b |.
(5.49)

As the friction velocity u τ increases faster with the Reynolds number than the conductive heat flux, the friction temperature decreases with the Reynolds number and so is the maximum temperature rms. The relative variation of temperature is

T max rms T (y + BL ) = T +,max rms (Re, Pr) sign(T b -T w )T + (y + BL (Re, Pr)) + Tw Tτ .
(5.50)

Introducing Eq. 5.49, it is found that

T max rms T (y + BL ) = T +,max rms (Re, Pr) sign(T b /T w -1)T + (y + BL (Re, Pr)) + α 0.5 β -1 Re 0.1 Pr 0.5 |T b /T w -1| -1
.

(5.51)

Finally, in air, the following result is obtained for the relative variation of temperature: where T + (y + BL = 12) ≈ 6.9 is calculated from the formula proposed by [START_REF] Kader | Temperature and concentration profiles in fully turbulent boundary layers[END_REF]. Equation 5.52 is plotted in Fig. 5.28 for several values of the ratio T b /T w which is less (resp. greater) than unity in heating (resp. cooling) systems. The level of temperature fluctuations decreases with the Reynolds number in all cases so that it can be expected to remain below 30% in all turbulent boundary layers in air. An additional section (Sec. 6.9) then presents results of the developed wall model combined with the turbulent Prandtl model proposed in chapter 5. Finally, as in direct numerical simulations, the bulk Reynolds number in LES is controlled by the method detailed in Appendix D.

T max rms T (y + BL ) = 2.5 6.9 sign(T b /T w -1) + 5.8 Re 0.1 |T b /T w -1| -1 , (5.52)

Introduction

Fully resolved Large Eddy Simulation (LES) of the inner layer of a wall-bounded turbulent flow requires highly resolved grids since the integral length scale becomes of the same order of magnitude as viscous scales in the close vicinity of the wall. The computational cost is then proportional to Re 2.4 [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF]. Hence, fully resolved LES is impracticable for wall-bounded flows at high Reynolds number, encountered in most of engineering applications, due to the prohibitive cost. Several kinds of approaches are commonly used in order to alleviate these difficulties: A wall model prescribes the correct wall shear stress to the LES that is too poorly resolved close to the wall to estimate it accurately. In hybrid RANS/LES, the simulation is switched from RANS in the inner layer to LES in the outer layer by the modification of the length scales [START_REF] Spalart | Comments of feasibility of les for wings, and on a hybrid RANS/LES approach[END_REF][START_REF] Nikitin | An approach to wall modeling in large-eddy simulations[END_REF] or the use of a blending function [START_REF] Hamba | A Hybrid RANS/LES Simulation of Turbulent Channel Flow[END_REF][START_REF] Shur | A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities[END_REF] for LES, the wall stress is estimated by using an algebraic wall function or by locally solving a simplified RANS equation. These approaches correspond to equilibrium-stress model and Two-Layer Model (TLM), respectively. The equilibrium-stress model has been firstly proposed by [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large reynolds numbers[END_REF]. It has then been successfully applied to turbulent channel flows and annuli flows [START_REF] Schumann | Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli[END_REF]. However, this model is restricted to simple flows, since it implies the existence of a logarithmic layer. In order to widen its use to more complex flows, the original equilibrium model has been modified by considering inclination of the elongated structure in near wall region [START_REF] Piomelli | New approximate boundary conditions for large eddy simulations of wall-bounded flows[END_REF]), pressure gradient [START_REF] Wang | LES with wall models for trailing edge aeroacoustics[END_REF][START_REF] Duprat | A wall-layer model for large-eddy simulations of turbulent flows with/out pressure gradient[END_REF], buoyancy [START_REF] Moeng | A large-eddy-simulation model for the study of planetary boundary-layer turbulence[END_REF] or chemistry [START_REF] Cabrit | Direct simulations for wall modeling of multicomponent reacting compressible turbulent flows[END_REF]. In two-layer models, turbulent boundary layer equations are resolved on a local embedded grid [START_REF] Balaras | Two-layer approximate References boundary conditions for large-eddy simulations[END_REF]. This approach has been extensively applied and assessed in different configurations [START_REF] Cabot | Approximate wall boundary conditions in the large-eddy simulation of high reynolds number flow[END_REF][START_REF] Wang | Dynamic wall modeling for large-eddy simulation of complex turbulent flows[END_REF][START_REF] Gungor | Direct simulation of subgrid turbulence in a high-Re wall-bounded flow[END_REF]. Moreover, since all wall modeled LES unavoidably suffers from the numerical and sub-grid error at the first grid point close to a wall [START_REF] Cabot | Approximate wall boundary conditions in the large-eddy simulation of high reynolds number flow[END_REF][START_REF] Nicoud | Large eddy simulation wall-modeling based on suboptimal control theory and linear stochastic estimation[END_REF]), an effective strategy has recently been proposed by [START_REF] Kawai | Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy[END_REF] to increase the accuracy of the information transmitted from LES in the outer layer to the inner thin turbulent boundary layer equations. A more detailed description of wall models for velocity is given in Refs. [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF][START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF][START_REF] Cabot | Approximate wall boundary conditions in the large-eddy simulation of high reynolds number flow[END_REF]. In order to deal with turbulent heat transfer and predict wall heat fluxes accurately, these wall models have to be extended to describe the thermal boundary layer as in [START_REF] Benarafa | Temperature wall modelling for large-eddy simulation in a heated turbulent plane channel flow[END_REF][START_REF] Rani | Boundary-Layer Equation-Based Wall Model for Large-Eddy Simulation of Turbulent Flows with Wall Heat Transfer[END_REF]). To the best of our knowledge, no wall model for LES has accounted for radiation effects, although radiation strongly modifies the temperature field in many applications, particularly in combustion processes at high pressure [START_REF] Gupta | Large-Eddy Simulation of Turbulence-Radiation Interactions in a Turbulent Planar Channel Flow[END_REF][START_REF] Soufiani | Radiation effects on turbulent heat transfer in channel flows of infrared active gases[END_REF][START_REF] Ghosh | Effects of radiative heat transfer on the structure of turbulent supersonic channel flow[END_REF]. It has been recently shown [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]) in coupled DNS-Monte Carlo simulations that radiation can significantly influence the temperature wall-law and the corresponding wall conductive heat flux. The temperature law is very different from the usual logarithmic law for strong radiation effects and has been observed to differ significantly under different radiative conditions. It is therefore unrealistic to hope for a general algebraic wall-law to account for these effects and a two-layer approach is then chosen. Besides, in order to predict the radiative field outside of the inner boundary layer, a reciprocal Monte-Carlo method is considered. The method is accurate and can be applied to complex geometries so that the proposed wall-model and its coupling with LES and the Monte-Carlo method remain general. The objective of this study is to account for radiation effects in the inner layer wall model to accurately predict wall stress and heat flux. Here, a two-layer model is retained where, in the outer layer, LES is coupled to a radiation Monte Carlo method as in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]. Coupled DNS-Monte Carlo results of Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF] LES wall model. The fluid and radiation models in both layers are detailed in Sec. 6.3, followed by a description of coupling between the inner and outer layers . Separate validations of the different model components are presented in Sec. 6.4. Finally, in Sec. 6.5, fully coupled results assess the model accuracy.

Wall-modeled LES coupled to radiation

In all fluid simulations, LES is here carried out in the outer layer and the boundary inner layer is modeled by solving 1D balance equations. For radiation, a reciprocal Monte Carlo approach is implemented to estimate the radiative power at all LES grid points and an analytical radiative 1D model is developed for the inner layer. For both radiation and fluid models, a particular care is brought to the boundary conditions, especially between the inner and outer layers. As explained in Ref. [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF], filtered equations in the inner layer are similar to averaged Navier-Stokes equations.Then, treating the unresolved inner layer [0, y w ] as a thin equilibrium boundary layer [START_REF] Wang | Dynamic wall modeling for large-eddy simulation of complex turbulent flows[END_REF][START_REF] Cabot | Approximate wall boundary conditions in the large-eddy simulation of high reynolds number flow[END_REF] leads to the following equations

Fluid model

d dy ( µ + µ t ) d{u || } dy = 0; d dy c p µ Pr + µ t Pr t d{T } dy + P R = 0, (6.1)
where the angled brackets • and curly brackets {•} denote Reynolds averaged and Favre averaged values respectively. u || is the tangential velocity, parallel to the wall, and y is the distance to the wall. T is the temperature, P R is the radiative power per unit volume. As in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF], µ, the dynamic viscosity, and c p , the thermal capacity at constant pressure, are functions of temperature, while the molecular Prandtl number Pr is set to 0.71. µ t and Pr t are the turbulent viscosity and the turbulent Prandtl number, respectively. µ t is computed using a mixing-length model [START_REF] Cabot | Near-wall models in large-eddy simulations of flow behind a backward-facing step[END_REF][START_REF] Cabot | Approximate wall boundary conditions in the large-eddy simulation of high reynolds number flow[END_REF])

µ t = κ ρ y τ w / ρ (1 -exp(-y + /A + )) 2 (6.2)
where κ = 0.4 and A + = 17 are here chosen in order to retrieve the velocity logarithmic law in the low-Reynolds DNS cases that are considered for validation. ρ is the gas density and τ w the shear stress at the wall. The wall coordinate y + is defined as

y + = ρ w u τ y µ w , u τ = τ w ρ w , (6.3)
where gas properties at the wall, ρ w and µ w , are used. In order to take into account the non-uniformity of the gas properties, an alternative semi-local coordinate y * (Huang et al. 1995) can also be used in Eq. ( 6.2) instead of y + , where local gas properties are considered

y * = ρ u * τ y µ , u * τ = τ w ρ , (6.4)
The mixing-length models based on wall coordinate y + and semi-local coordinate y * will subsequently be referred as standard and semi-local mixing-length models, respectively. Three models for the turbulent Prandtl number are considered in the following. First, a constant value Pr t = 0.9 which is a usual crude assumption in wall models. Secondly, in order to account for variation of Pr t in the boundary -5.165 µt / µ ))

. (6.6)

Outer layer LES

A low-Mach code, YALES2 (Moureau et al. 2011c;Moureau et al. 2011b), is used for the Large Eddy Simulation (LES) of the main flow. The same numerical set-up as in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]) is retained: A 4th-order central difference finite-volume scheme with a 4th-order time integration scheme.

Under the low-Mach number approximation, the spatially filtered instantaneous mass, momentum and energy balance equations write

∂ρ ∂t + ∂(ρ u i ) ∂x i = 0, (6.7)
∂(ρ u i ) ∂t + ∂(ρ u i u j ) ∂x j = - ∂τ SGS ij ∂x j - ∂p ∂x i + ∂τ ij ∂x j + S i , (6.8)
∂(ρ h) ∂t + ∂(ρ u j h) ∂x j = - ∂q SGS j ∂x j + ∂p ∂t - ∂q cd j ∂x j + P R , (6.9)
where p is the pressure, u i the velocity components and • and • denote filtered and mass-weighted filtered quantities in the context of LES. The enthalpy per unit mass h is expressed as h = ∆h 0 + T T 0 c p (T ′ ) dT ′ , where c p is the mixture thermal capacity at constant pressure, T 0 a reference temperature and ∆h 0 the corresponding standard formation enthalpy. τ ij and q cd i are the viscous shear stress tensor and the conductive heat flux vector respectively. S i is a driving force source term to obtain the intended bulk Reynolds number in channel flow simulations. Based on the Sub-Grid Scale (SGS) eddy-viscosity concept, the SGS stress tensor τ SGS ij = -ρ( u i u j -u i u j ) and the SGS heat flux q SGS j = -ρ( hu j -h u j ) are modeled as where S ij is the filtered shear stress tensor. The Sigma model [START_REF] Nicoud | Using singular values to build a subgrid-scale model for large eddy simulations[END_REF]) is used for modeling the SGS kinematic viscosity ν SGS . The SGS thermal conductivity λ SGS is computed from the SGS Prandtl number Pr SGS which is set to 0.9 λ SGS = ρ c p ν SGS Pr SGS . (6.12)

τ SGS ij - 1 3 δ ij τ SGS kk = -2ρν SGS (S ij - 1 3 δ ij S kk ), (6.10)
q SGS j = -λ SGS ∂T ∂x j , (6.11)

Radiation model

For the radiation field, a Monte-Carlo method is used to calculate the radiative power field in the outer layer, i.e. at all the LES grid points. For the inner layer, the radiation field is analytically obtained from a one-dimensional model that uses the intensity field obtained from the Monte Carlo approach as a boundary condition for the point with wall distance y w , where information from the LES grid is fed back to the wall model.

Monte Carlo approach

The general organization of the radiation model, based on a reciprocal Monte Carlo approach, has been detailed by [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF]. The precise approach used here is the Optimized Emission-based Reciprocity Monte Carlo Method (OERM) [START_REF] Zhang | Optimized emission-based reciprocity monte carlo method to speed up computation in complex systems[END_REF], as in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]. This method allows the convergence to be locally controlled while it overcomes the drawback of the original Emission-based Reciprocity Method [START_REF] Tessé | Radiative transfer in real gases using reciprocal and forward monte carlo methods and a correlated-k approach[END_REF] in the cold region and greatly increases the computational efficiency. This method also allows the spectral anisotropic radiation intensity field I - ν (y w , µ) at any point of abscissa y w of the LES grid, required by the inner layer radiation model, to be determined, as detailed in Sec. 6.3.2.2 and Appendix. The one-dimensional configuration for radiation in the inner layer is shown in Fig. 6.2. The opaque wall is isothermal at temperature T w and its emissivity ǫ ν

I - ν (y, µ 2 ) θ 1 θ 2 T w ǫ ν
I - ν (y, µ)2πdµ + 1 0 I + ν (y, µ)2πdµ dν (6.15)
where ν is the radiation wave number, µ is the cosine of the angle θ defined in Fig. 6.2, and I • ν (y) is the equilibrium spectral intensity at the temperature T (y) associated with the current position y. The current spectral intensity in a positive y direction I + ν (y, µ), associated with µ > 0, and the spectral intensity in a negative direction I - ν (y, µ), associated with µ < 0 are given by

I + ν (y, µ) = τ ν0y (µ)I + ν (0) + y 0 κ ν (y ′ )I • ν (y ′ )τ νy ′ y (µ)
dy ′ µ with: y ′ < y; dy ′ > 0; µ > 0, (6.16)

I - ν (y, µ) = τ νywy (µ)I - ν (y w , µ) + y yw κ ν (y ′ )I • ν (y ′ )τ νy ′ y (µ)
dy ′ µ with: y ′ > y; dy ′ < 0; µ < 0, (6.17) where (6.18) and τ ′ νy ′ y is the directional spectral transmissivity between y ′ and y given by Here e ν (y ′ , y)/µ > 0 is the spectral optical thickness between y ′ and y in the direction µ while κ ν is the spectral absorption coefficient of the medium. Note that: i) dy ′ and µ in Eqs. (6.16) and (6.17), and e ν (y ′ , y) and µ in Eq. (6.19) are simultaneously positive or negative; ii) As the wall reflection is assumed diffuse, the intensity leaving a wall I + ν (0) is isotropic, whereas I - ν (y w , µ) the intensity entering the inner layer at abscissa y w is anisotropic. By introducing the exponential integral function, detailed in Ref. [START_REF] Modest | radiative heat transfer[END_REF], i.e.

I + ν (0) = ǫ ν I • ν (T w ) -2 (1 -ǫ ν ) 0 -1 I - ν (0, µ) µdµ,
E n (X) = In this analytical radiation model in the boundary inner layer, the radiative power P R (x) is a function of the temperature field within the fluid inner layer model given by Eq. (6.1). It is then worth noticing that the effects of turbulence fluctuations on the radiative power are not accounted for within the inner layer only. In fact, these effects can indeed be neglected in channel flows with nonreacting gases where fluctuations of temperature remain moderate, as shown by Ref. [START_REF] Gupta | Large-Eddy Simulation of Turbulence-Radiation Interactions in a Turbulent Planar Channel Flow[END_REF]) and by post-processing of DNS results in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF] .

Coupling of inner and outer layers radiation and turbulence models

The purpose of wall-modeled LES is to overcome the under-resolution of the boundary inner layer, which leads to erroneous estimations of wall temperature and velocity gradients. For each LES grid point on the wall and each time step, the set of equations (6.1) combined with the analytical radiation model is solved with an iterative procedure and provides an accurate estimation of the wall stress τ w and wall conductive heat flux q cd w . The procedure used to For each grid point on the wall, LES provides velocity u yw and temperature T yw at a distance y w from the wall to the 1D turbulence model. It also provides the resolved temperature field T LES to the Monte-Carlo method.

The Monte-Carlo method calculates both the radiative power P R MC , source term of the energy balance equation in LES, and the incoming intensities I ν (y w , [µ j , µ j+1 ]) required by the 1D radiation model. Note that in Monte Carlo model, when a ray enters the inner layer, the 1D model temperature field is used to determine the exchanged energy between the point initiating the ray and the inner layer. The accuracy of the results is then improved when compared to results associated with the LES temperature field (much coarser in the near wall region).

LES

Monte-Carlo

T LES P R MC τ w , q cd w u yw , T yw T 1D I - ν (y w , [µ j , µ j+1 ]) 1D Wall Model Fluid Radiation P R 1D T 1D Figure 6
.3: Scheme of coupling between LES, the Monte-Carlo method and wall model. T 1D and P R 1D are the inner layer temperature and radiative power fields on the embedded grids.

Separate validation of the turbulence and radiation wall model components

The previously described 1D turbulence wall model is a priori validated against DNS data from Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF] accounting or not for radiation. This a priori study consists in comparing 1D profiles of average quantities predicted by the wall-model while taking some information from DNS cases. DNS cases from Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]) that are retained is this study are defined in walls T w,h and T w,c are different, as shown in Fig. 6.4, and radiation is not accounting for. For cases A, B and C, the corresponding cases with radiation are named A_RAD, B_RAD and C_RAD respectively and opaque wall emissivities are set to 0.8. Two different cases at one atmosphere, named D_RAD03 and D_RAD08, are considered for case D with wall emissivities of 0.3 and 0.8 respectively. The gas radiative properties are modeled by using the CK model [START_REF] Soufiani | High temperature gas radiative property parameters of statistical narrow-band model for H 2 O, CO 2 and CO, and correlated-K model for H 2 O and CO 2[END_REF] for atmospheric pressure cases and its weak absorption limit for cases at 40 atm. The domain is defined in Fig. 6.4 with : δ = 0.1m, L X = 2πδ, L Y = 2δ, L Z = πδ. A non-reacting CO 2 -H 2 O-N 2 gas mixture is considered and the molar fractions of CO 2 , H 2 O and N 2 are set to 0.116, 0.155 and 0.729 to mimic combustion exhaust gases. Similarly to Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF], the dynamic viscosity µ and thermal capacity c p are computed as functions of temperature from the CHEMKIN package [START_REF] Kee | A fortran computer code package for the evaluation of gas-phase, multicomponent transport properties[END_REF][START_REF] Kee | CHEMKIN-II: A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase Chemical Kinetics[END_REF] for the chosen mixture composition. The thermal conductivity λ is computed from the Prandtl number Pr= 0.71.

In the following sections, different wall-model components are separately validated.

Validation of the turbulence wall model without radiation

An integration of Eq. (6.1), in which the wall friction stress and wall heat flux results of DNS are imposed as boundary conditions, leads to velocity and temperature profiles in the near wall region. Since the proposed mixing-length model is only valid in the inner layer, only results within the near wall region (0, 0.2δ) are presented. Velocity profile of case B that presents the most important wall temperature and gas properties variations is shown in Fig. 6.5. predicted by the wall model based on standard mixing-length model, Eq. ( 6.3). The semi-local model, Eq. (6.4), enables to correct this effect. When temperature variations are negligible, both models are identical and perform accurately as shown in Fig. 6.5 for the DNS case from Ref. [START_REF] Hoyas | Scaling of the velocity fluctuations in turbulent channels up to Re τ = 2003[END_REF] characterized by constant temperature and a larger Reynolds number. Regarding the temperature profile, the two mixing-length models and the three models of turbulent Prandtl number defined in Sec. 6.3.1.1 are investigated in Fig. 6.6. Best results are obtained using the semi-local mixing-length model and Eq. (6.6) for the turbulent Prandtl number. These modeling conditions are henceforth retained. Velocity and temperature profiles for cases A and C are shown in Fig. 6.7. Very good agreement with DNS data is observed.

Validation of the turbulence wall model with prescribed radiative power field

In this validation cases, the reference radiative power issued from the radiation Monte Carlo method coupled to DNS is also used in addition to the DNS wall friction stress and wall heat flux in Eq. (6.1). Figure 6.8 shows that the complete wall model (semi-local mixing-length model, Eq. (6.6) for Pr t and radiative power source term) accurately predicts the temperature field for the two near wall regions of the three cases A_RAD, B_RAD and C_RAD. Note that the results of wall models that do not account for radiative power source term strongly deviate from the DNS corresponding results, which indicates that in these three cases, radiation strongly modifies the mean temperature field within the inner layer.

Validation of the radiation wall model with prescribed temperature field

In order to validate the radiative analytical wall model, Eq. (6.22) is solved using the mean temperature profile obtained from DNS data in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF] and N µ spectral anisotropic incoming intensity values issued from the reference cases post-processing as explained in Appendix. The influence of the number of angular sectors N µ and of the wall distance y w where anisotropic intensities are calculated is investigated. Results are shown in Fig. 6.9. The radiative power is accurately predicted with N µ = 2 at each value of y w between 0.05δ and 0.2δ. N µ = 2 is therefore used in the following. Note that results 6.5 Results of the turbulence and radiation wall models fully coupled with LES and Monte Carlo method

In this section, the turbulence and radiation wall models are fully coupled to LES and Monte Carlo models for the same cases A to D as in the previous section without prior knowledge from the reference DNS cases. For all these cases, the LES grid consists of 36 × 36 × 36 points and the first off-wall point locates at y + = 15, while all the other grid spacings are uniform. The wall model, defined by Eqs. (6.1, 6.22), is resolved at each point on the wall on a local embedded grid (see Fig. 6.1). Given first guesses of the wall shear stress and wall heat flux, this set of equations is iteratively solved while ensuring grid convergence.

Cases without radiation

When coupling LES and wall models, numerical and sub-grid model errors at the first few points near the wall induce perturbations within the LES domain [START_REF] Cabot | Approximate wall boundary conditions in the large-eddy simulation of high reynolds number flow[END_REF][START_REF] Nicoud | Large eddy simulation wall-modeling based on suboptimal control theory and linear stochastic estimation[END_REF], leading to the so-called "logarithmic law mismatch" [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF]. The term "logarithmic law mismatch" specifically refers to cases where a logarithmic law is found for the velocity profile. A similar mismatch between the wall modeled LES and the reference results can also be observed for the temperature profile of the different cases studied here, where the temperature law is not logarithmic as soon as radiation is accounted for.

In order to reduce this error, [START_REF] Kawai | Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy[END_REF] have proposed to place the feed back point of the wall model, placed at the wall distance y w , further away from the wall than the first off-wall LES grid point. Following this approach, different y w positions are compared, i.e. y w located at the first point, second and third off-wall points, denoted case Y1, Y2 and Y3 respectively. The corresponding embedded grids, which are stretched along the wall normal direction, contain 30, 60 and 100 points. Fig. 6.10 shows that, for case A, the case Y3 presents the best agreement for both u + and T + profiles as expected. This trend is retrieved for wall friction stress τ w and conductive heat flux q cd w values shown in Tab. 6.2. The feed back point position does not affect the root mean square (rms) profiles of u, v, w and T (see Fig. 6.11). Henceforward, the third off-wall point is retained as a feed back point. Results for cases B and C are presented in Fig. 6.12 and Tab. 6.2 demonstrating the relevance of this choice. Table 6.2: Comparison of mean friction stress τ w and conductive heat flux q cd w for cases A, B and C between wall-modeled LES (Y1, Y2, Y3) and DNS results from Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]. For case A, Y1, Y2 and Y3 correspond to wall-normal distances of 0.05δ, 0.098δ and 0.156δ respectively.

Cases involving radiation at high pressure

The finally proposed wall model based on semi-local mixing-length model, Eq. (6.6) for Pr t , the radiation analytical model described in Sec. 6.3.2.2 and using the third off-wall point for coupling (case Y3) will be denoted as "new" wall model. In order to demonstrate the accuracy of this model, it is compared to a "standard" wall model based on standard mixing-length model and constant turbulent Prandtl number, that does not account for radiation effects within the inner layer. Note that this "standard" model is also coupled at the third off-wall point and that radiation is treated with the Monte-Carlo method in the outer layer. Results without any wall model are also presented. DNS results of u + and T + in cases A_RAD, B_RAD and C_RAD are compared in Fig. 6.13 to those predicted by wall-modeled LES with the "standard" and "new" approaches. The velocity field is similarly reproduced by both models whereas the "new" model shows a significant improvement of the predicted temperature profile. Values of τ w and q cd w for the different cases are given in Tab. 6.3 where the same conclusion is retrieved.

Case involving radiation at 1 atm

Cases A to C correspond to optically thick media at 40 atm, in which the effects of radiation are important. Case D, considered in this paragraph and defined in Tab. 6.1, deals with media at atmospheric pressure. LES of the two cases D_RAD03 and D_RAD08, with wall emissivities set to 0.3 and 0.8 respectively, with the "new" and "standard" wall models are compared with DNS data from Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]). Figure 6.14.a compares the u + profiles that are only shown on the cold side (similar on the hot side) and demonstrates that both models predict the correct result as in the high pressure cases. The same behavior is observed for the wall shear stress in Tab. 6.3. Regarding the T + profiles (Figs. 6.14.b and 6.14.c) and wall conductive fluxes given in Tab. 6.3, the "new" model leads to the best agreement with DNS-Monte Carlo results compared to the "standard" model. However, since the radiation effects are weaker at 1 atm than at 40 atm, the difference between the results of the two wall-modeled LES are smaller than those at high pressure, especially on the cold side.

Conclusion

A wall model for LES has been proposed to take into account the effects of radiation within the turbulent boundary inner layer. The model follows a two-layer approach with an embedded grid for each point at the wall. It is composed of 1D equilibrium thin boundary layer equations that are solved on each embedded grid along with mixing-length and turbulent Prandtl models and an analytical expression of the radiative source term within the inner layer. First results have shown that: i) A semi-local scaling of the mixing-length is necessary to account for the variations of gas properties, such as mass density and dynamic viscosity, within the inner layer; ii) A non-constant turbulent Prandtl number, here Eq. (6.6), gives better results; iii) In the studied cases, the mean incoming intensity at the outer boundary of the wall model has required the discretization of the corresponding solid angle into at least two angular sectors to accurately account for anisotropy; iv) The point where information is transmitted from LES to the wall model has been moved away from the wall while remaining in the validity range of the wall-model. Here, the third grid point off the wall has been used in order to decrease the mismatch in the obtained wall laws due to numerical and sub-grid errors in the wall vicinity. The LES is coupled to a reciprocal Monte-Carlo method as it could be done in any configuration to tackle radiation and turbulent flows. Here, the combination of the proposed wall model and LES/Monte-Carlo method is compared to DNS/Monte-Carlo results on turbulent channel flows. Coupled results are compared to those obtained with a standard wall model which does not account for radiation effects within the inner layer. Note that this standard model with LES and Monte-Carlo simulations is already quite advanced to accurately predict wall heat transfer in the presence of radiation. Nonetheless, coupled results for high pressure cases show that the new wall-model greatly improves the accuracy of the predicted results compared to a standard wall model. The same trend is retrieved in atmospheric cases even though the magnitude of difference is smaller than the one at high pressure. It has been finally shown that accounting for radiation effects within the modeled boundary inner layer is necessary in several conditions, especially in high pressure applications, as soon as the radiative energy transfer is strong enough to modify the temperature wall law. Table 6.3: Comparison of mean friction stress τ w , wall conductive heat flux q cd w and total wall radiative flux q R w of cases with radiation between wall-modeled LES (new, standard and no wall model) and DNS results from Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]). ), boundary condition of the wall model, is computed from the reciprocal Monte-Carlo method by using a virtual wall. µ j and µ j+1 correspond to the cosine of angles θ j and θ j+1 respectively.

As described in Sec. 6.3.2.2, the proposed wall model requires mean anisotropic spectral incoming intensities I ν (y w , [µ j , µ j+1 ]) at a grid point y w averaged over the solid angles corresponding to the ranges [µ j , µ j+1 ] as a boundary condition. The emission-based reciprocity Monte Carlo method defined in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]) is used in this paper. Consequently, the mean spectral intensities can easily be deduced from reciprocal emission phenomena issued from all the other cells of the system. The precise used procedure is the following: i) The grid point at wall distance y w is assumed to be a virtual opaque wall of temperature T (y w ) and of arbitrary wall emissivity ǫ = 1; ii) For each range [µ j , µ j+1 ], a large number of shots are emitted in the Monte Carlo method from this virtual wall in a randomly determined direction with the cosine of polar angle within the range [µ j , µ j+1 ], and in any elementary spectral range ∆ν ; iii) The reciprocity method allows Φ ∆ν (y w , [µ j , µ j+1 ]), the contribution to the radiative flux at the virtual wall associated with both [µ j , µ j+1 ] and ∆ν, to be determined; vi) The absorbed flux at the wall Φ a ∆ν (y w , [µ j , µ j+1 ]), and consequently the mean spectral absorbed intensity I a ν (y w , [µ j , µ j+1 ]), associated with [µ j , µ j+1 ] and ∆ν is then calculated as

Φ a ∆ν (y w , [µ j , µ j+1 ]) = Φ ∆ν (y w , [µ j , µ j+1 ]) + Φ e ∆ν (y w , [µ j , µ j+1 ]), (6.23) I a ν (y w , [µ j , µ j+1 ]) = Φ a ∆ν (y w , [µ j , µ j+1 ]) 2π µ j+1 µ j µ dµ ∆ν = Φ ∆ν (y w , [µ j , µ j+1 ]) 2π µ j+1 µ j µ dµ ∆ν +I 0 ν (y w ); (6.24)
v) The mean spectral incoming intensity I ν (y w , [µ j , µ j+1 ]) required by the wall model is then equal to I a ν (y w , [µ j , µ j+1 ]), as the virtual wall absorptivity has been set to unity. Only the most challenging case B_R is retained in this test since it reveals the most important radiation effects and gas property variations. Following the methodology of the wall model validation, a similar a priori validation taking as little knowledge from DNS as possible. The test has the following features: i) Only the energy equation of Eq. (6.1) is solved and the mean temperature at position y w ≈ 0.2δ is taken from DNS and used as boundary condition; ii) The turbulent viscosity µ t is computed from the mixing-length model (Eq. ( 6.2)) along with the semi-local coordinate y * (Eq. ( 6.3)), using the DNS results for the wall shear stress τ w ; iii) An iterative procedure is carried out to solve the energy equation; iv) The radiative power term in Eq. (6.1) is interpolated from DNS results and is fixed during the iterative procedure; v) The turbulent Prandtl number is calculated either from Eq. (6.6) or from the new model Eq. (5.45).

The results of the wall model are compared with DNS datas from [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF] in Fig. 6.16. On the cold side of case B_R, except in the very vicinity of walls where turbulent heat transfer is negligible compared to conduction, the turbulent Prandtl number profile predicted by the new model Eq. (5.45) agrees well with DNS datas whereas it is underestimated by Eq. (6.6). However, profiles of mean temperature T + show only small differences between results obtained with both turbulent Prandtl number models. In fact, in this case as for the other ones that have been studied, the radiative transfer dominates over the turbulent heat transfer and hence the accuracy of turbulent Prandtl number is not critical. However, when turbulent convection is predominant or is of the same order of magnitude as radiative transfer, the influence of turbulent Prandtl number model is expected to be more important, although radiation effects might then be weaker.

Regarding the Pr t profile on the hot side, the accuracy of the predicted results is also improved by the new model. However, a noticeable gap is still observed between wall model results and DNS data, indicating the necessity of a further improvement. Nevertheless, as for the cold side, the profile of T The obtained criterion results are presented and detailedly analyzed to investigate a large set of conditions, which would not be feasible with DNS or LES.

Description of the criteria

Since calculating the radiation field is generally very expensive, it could be useful if one could predict the importance of radiation before performing an accurate simulation. In this section, several criteria are proposed in order to determine whether it is necessary to simulate the radiation field and further, to predict if the wall model proposed in chapter 6 is required under the studied flow conditions. pared to the conductive one, it is reasonable to define a criterion as following:

Criteria for radiative flux

R q,coupled = q R * w ({T coupled }) q cd w ({T coupled }) (7.1)
where subscript coupled denotes quantities resulting from a simulation performed with a full coupling between fluid mechanics and radiation while {•} relates to the whole field of quantities. q R * w ({T coupled }) and q cd w ({T coupled }) are the wall radiative and conductive flux respectively. Like that in section 4.3.5 of chapter 4, the flux exchanged between the walls are not included in the the wall radiative flux q R * w considered here. However, since the objective of the criterion is to predict the importance of radiative flux before performing a real coupled radiation/flow field simulation, the quantities q R * w ({T coupled }) and q cd w ({T coupled }) are generally not available. To resolve this problem, the mean temperature field T obtained without accounting for radiation is then used to compute a fictitious wall radiative flux q R * w ({T }). And another criterion can then be defined with this radiative flux q R * w ({T }) and the corresponding wall conductive heat flux q cd w ({T }) as following:

R q = q R * w ({T }) q cd w ({T }) . (7.2)
These two criteria R q and R q,coupled have been compared for all the channel flow cases presented in chapter 4 and both walls (cold and hot). These results are presented in Fig. 7.1 (a) It reveals that for all the cases, values of the two criteria are highly correlated, indicating that R q can be used to roughly predict the contribution of the wall radiative flux to the total wall heat flux instead of R q,coupled .

Criteria for wall-scaled temperature

Since radiation can strongly influence the temperature law of turbulent boundary layer under certain conditions (as presented in chapter 4 and 5), a wall model considering this effect has been proposed and validated in chapter 6. In this part, a criterion is proposed to predict if it is mandatory to use the new wall model or not in order to have a good estimation of the temperature field in the boundary layer.

In near wall regions, when omitting viscous heating and unsteadiness, the mean balance equation of energy then writes

d dy c p coupled µ coupled Pr + µ R t,coupled Pr R t,coupled dT coupled dy +P R ({T coupled }, y) = 0 (7.3)
where subscript coupled also denotes quantities associated with the case with radiation. P R is the radiative power per unit volume, µ and c p are the dynamic

dy + = 1 - ∆q R ({T coupled }, y) q cd w ({T coupled }) , (7.4)
where q R ({T coupled }, y) is the radiative flux and

∆q R ({T coupled }, y) = - y 0 P R ({T coupled }, y ′ )dy ′ (7.5) = q R ({T coupled }, y) -q R w ({T coupled }) (7.6
) is the change of radiative flux. The non-dimensional turbulent viscosity is defined as: µ + t = µ t /µ. The non-dimensional distance y + , based on the distance to a wall y and the non-dimensional temperature T + write

y + = ρ w y u τ µ w , T + = |T -T w | T τ , (7.7)
where T τ is the friction temperature which is defined by

T τ = |q cd w | ρ w c pw u τ , u τ = τ w ρ w 1/2 . (7.8)
and τ w is the wall shear stress.

For cases where no radiation is considered, Eq. (7.4) turns into

1 Pr + µ + t Pr t dT + dy + = 1. (7.9)
Assuming that the change of µ + t /Pr t due to radiation is not important, the difference of wall-scaled temperature between the case without/with radiation from Eq. (7.4) and Eq. (7.9) can be approximated as

1 Pr + µ + t Pr t d(T + -T + coupled ) dy + ≈ ∆q R ({T coupled }, y) q cd w ({T coupled }) (7.10)
Note that the quantities without a subscript coupled relates to those of cases without radiation. An integration of Eq. (7.10) leads to is defined as

T + -T + coupled T + ≈ 1 T + y + 0 1 1 Pr + µ + t Prt ∆q R ({T coupled }, y) q cd w ({T coupled }) dy ′+ . ( 7
∆ ref t + ,coupled y=0.1δ =   1 T + y + 0 1 1 Pr + µ + t Prt ∆q R ({T coupled }, y) q cd w ({T coupled }) dy ′+   y=0.1δ
.

(7.12)

The wall-normal distance 0.1δ is approximatively the limit of the boundary inner layer, where the upcoming turbulent models are valid. The values of this criterion

∆ ref t + ,coupled y=0.1δ
for all the cases from chapter 4 (except the case with different emissivities on the two walls) are compared against DNS results of

T + -T + coupled T + y=0.1δ
in Fig. 7.1 (b). It shows that criterion values calculated from Eq. (7.12) are nearly the same as DNS results, despite the assumption that µ + t /Pr t is not modified by radiation, which is not entirely true for the turbulent Prandtl number as shown in chapter 6. Similarly, since results of cases with radiation (i.e. ∆q R ({T coupled }, y), q cd w ({T coupled }) etc.) are generally not available, another criterion based only on the flow field of the case without radiation is defined as (7.13) where q cd w ({T }) is employed instead of q cd w ({T coupled }) and radiative flux q R ({T }, y) is calculated from the mean temperature field of the case without radiation. The results of this criterion are also shown in Fig. 7.1 (b). Although the agreement between the value of this criterion and

∆ ref t + y=0.1δ =   1 T + y + 0 1 1 Pr + µ + t Prt ∆q R ({T }, y) q cd w ({T }) dy ′+   y=0.1δ ,
T + -T + coupled T + y=0.1δ
is poorer, the trend is correct. A higher criterion value implies a larger change of T + due to radiation. Hence, the criterion

∆ ref t + y=0.1δ
based only on the case without radiation can be used.

If the criterion R q tells us when a radiation simulation is needed, then ∆ ref

t + y=0.1δ
indicates whether the simulations of radiation and flow field need to be fully coupled, and more importantly, whether the wall model considering radiation is required to predict wall heat flux accurately. The values of these two criteria in different turbulent channel flows will be calculated in the following section.

Numerical model for flow field

These two criteria will later be used to perform a parametric study in channel flow in order to identify conditions where simulating the radiative transfer is important. To obtain the velocity field in these channel flows, the k-ǫ model is employed, where the turbulent viscosity is also obtained. An empirical formula of turbulent Prandtl number is used to resolve the mean temperature field.

Main equations

Firstly, to simulate the mean velocity field, the k-ǫ model from Ref. [START_REF] Abe | A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows-ii. thermal field calculations[END_REF] . The treatment of the source terms, S u and S T , along with special attentions dedicated to the discretization and stabilization of the main equations are given in Appendix E. Once the temperature field T is obtained, the radiative flux q R ({T }, y) can then be computed analytically [START_REF] Taine | Transferts thermiques : Introduction aux transferts d'énergie[END_REF] in the whole channel. This is similar to the radiation wall-model in chapter 6 which takes into account anisotropic boundary conditions, although, here, the classical textbook solution is used since radiative boundary conditions are located at the two channel walls and are then isotropic.

Validation of flow field

Using the set of equations described above, the channel flow case of Kim and Moin (Kim et al. 1987;Kim and Moin 1987) is computed and the results for mean velocity u and temperature T are compared in Fig. 7.2 (here the temperature acts as a passive scaler and only the results of the passive scalar with Pr=0.71 are shown). It reveals that both the velocity and temperature are predicted with a very good agreement with corresponding DNS results validating the numerical strategy.

Results of the criteria

The aforementioned numerical approach is used to simulate many channel flows for a large range of hydraulic diameters of channel, pressure, wall emissivity and practical Reynolds number. For all the cases under consideration, the temperatures at walls and in the center are set to typical values in gas turbine combustors, i.e., 800K and 1800K respectively. The medium is non-reacting CO 2 -H 2 O-N 2 gas mixture, as that used in previous studies of this thesis. Molar fractions of CO 2 , H 2 O and N 2 are 0.116, 0.155 and 0.729. The dynamic viscosity µ is computed as a function of temperature from the CHEMKIN package [START_REF] Kee | A fortran computer code package for the evaluation of gas-phase, multicomponent transport properties[END_REF][START_REF] Kee | CHEMKIN-II: A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase Chemical Kinetics[END_REF] for the chosen mixture composition. The thermal conductivity λ is computed from a Prandtl number Pr. In all considered configurations, the Prandtl number is very close to the chosen value 0.71. The obtained results of the two criteria are presented and analyzed in this section. . 7.3 (a) shows the contour plot of criterion R q when practical Reynolds number Re D h and hydraulic diameter of channel vary. The wall emissivity is 0.8 and the pressure is set to 40 atm in order to have relatively optically thick medium. a criterion value of 0.2 is chosen as a threshold to start a radiation calculation since a weight of 20% of radiative flux is already important. As expected in such a case, the region where R q is higher than 0. ≥ 0.2. Obviously, this gray region is smaller than the one of R q , indicating that in some region, although the wall radiative flux is important, the radiation effect is not strong enough to influence T + profiles near the wall.

Variation of radiation effect at high pressure

The magnitude of criterion value generally increases with the hydraulic diameter of channel while the change with Reynolds number is more complicated. In fact, regarding the influence of Reynolds number on wall-scaled temperature, two aspects have been identified in Ref. [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]): i) As the Reynolds number increases, the weight of turbulent convective heat transfer increases and hence it yields a relatively less important radiation effect. Consequently, T + profile approaches to the usual wall law, i.e., absolute value of

∆ ref t + y=0.1δ
decreases; ii) The temperature gradient in vicinity of walls increases with Reynolds number. Gas temperature at a given position then differs more from wall temperature while the difference with the average gas temperature is reduced. Consequently, the importance of gas-gas radiation declines while gaswall radiation effect is enhanced. Moreover, as indicated in Ref. (Zhang et value whereas gas-wall radiation has an opposite effect. Therefore, as Reynolds number increases, value of the criterion is expected to decrease and even become negative when gas-wall radiation dominates over gas-wall radiation. Indeed, in Fig. 7.3 (b), it turns out that when Reynolds number is relatively small (smaller than 40000), gas-gas radiation is dominant and a positive criterion value is obtained. Both of the two aforementioned aspects tend to decrease the criterion value when Reynolds number increases. Furthermore, if Reynolds number becomes even higher, a negative criterion value appears indicating that gas-wall radiation overcomes gas-gas radiation. However, this negative value is small since radiation effect is weak at high Reynolds number. Figure 7.4 shows the contour plot of the two criteria when practical Reynolds number and wall emissivity change for a pressure set at 40 atm and a hydraulic diameter of channel fixed at 0.4 m. It is observed as expected that the radiative flux, and hence R q , increases with wall emissivity since the gas-wall radiation is enhanced. On the other hand, as wall emissivity increases, the weight of gasgas radiation decreases and hence the value of criterion ∆ ref

t + y=0.1δ
decreases from a high positive value to even a negative value when with a high Reynolds number.

The contour plots of the criteria associated with variable hydraulic diameter and wall emissivity are shown in Fig. 7.5. The pressure condition is the same as before while practical Reynolds number is 48,000. As described above, both these two criteria increase with hydraulic diameter of channel while the trends with wall emissivity are opposite. Therefore, for criterion R q , the gray region is on the top-right corner while it is located on the bottom-right corner for

∆ ref t + y=0.1δ
.

Variation of radiation effect at 1 atm

The results of criteria at 1 atm are demonstrated in Fig. 7.6, 7.7 and 7.8 when changing different flow conditions. The trends of these results are similar to those at 40 atm whereas the magnitude of the values and the size of the gray region decrease owing to decrease in the optical thickness at 1 atm compared to that at 40 atm. Unsurprisingly, smaller radiation effects are then observed at 1 atm.

Variation of radiation effect with pressure

Since the gas absorption coefficients increase and hence the optical thickness of medium, with pressure condition, stronger radiation effects are achieved as the pressure increases. This is observed in Figs. 7.9, 7.10 and 7.11 where both criteria values are plotted. These observations under a multitude of conditions have confirmed and generalized the conclusions obtained from previously analyzed DNS data: Due to an increased optical thickness of medium, radiation effect is enhanced and higher criterion values are obtained as hydraulic diameter of channel or/and pressure increase. When wall emissivity changes, the trend of the two criteria are opposite since the wall radiative flux is mainly due to gas-wall radiation while change of wall-scaled temperature, on the other hand, is generally dominated by gas-gas radiation. The two criterion values generally decrease with Reynolds number, while a small negative value of the criterion for wall-scaled temperature can be achieved when gas-wall radiation overcomes gas-gas radiation in flows characterized by high Reynolds number. Finally, the proposed criteria and their evaluation method form a practical tool for engineers to determine whether radiative effects must be accounted for.

"thesis_yufang" -2013/11/27 -19:31 -page 189 -#189

✐ ✐ ✐ ✐ ✐ ✐

Conclusion

The general objective of this thesis is to assess a better understanding and prediction of the influence of radiation on turbulent flows. Although various studies have been dedicated to this problem, there was still no general trend or understanding of radiation effects in turbulent boundary layers. Therefore, a detailed physical study of such effects on the mean flow field and fluctuation field in turbulent channel flows has been performed in this thesis. A new wall model for LES has also been developed to account for radiation in the near wall region. Finally, two criteria have been proposed to assess the necessity of performing a radiation simulation and to predict whether the proposed wall model is required under a given flow condition. The main achievements of this thesis are described in the following.

Main achievements

• For radiation simulations in this work, the Monte-Carlo method is retained in order to obtain high-fidelity results. Among various approaches of Monte-Carlo methods, the Emission-based Reciprocity Method (ERM) is more efficient since it enables to control the uncertainty of the results locally or/and to shot photon bundles only from the points where results are required. However, this method is inaccurate in the treatment of absorption in cold regions. In order to overcome this drawback of ERM, an Optimized Emission-based Reciprocity Method (OERM) is proposed and validated in the first part of this thesis. OERM applies a new frequency distribution function based on the maximum temperature of the whole domain while it still fulfills the real emission distribution function of any temperature by adding a corrective factor to each shot. This method has been validated in several benchmark cases by a comparison with analytical reference solutions and results of some other reciprocity Monte-Carlo methods when available.

• In the second part of this thesis, the coupling between turbulent convection and radiation in channel flows has been investigated for different temperature, optical thickness (pressure) and wall emissivity conditions. Direct Numerical Simulation has been performed for the flow field and has been fully coupled to radiation simulation using OERM combined to the CK model to describe radiative properties. In the studied conditions, the structure of the temperature field and the wall conductive flux often strongly differ from results without radiation. Classical wall log-laws for temperature are then modified by the global radiation effects. It has been understood that gas-gas and gas-wall radiation interactions generate opposite effects on the temperature and flux fields. The first one tends to increase wall conductive flux while the second one to decrease it. Finally, the averaged temperature profiles and wall conductive fluxes were shown to strongly depend on: i) pressure, that controls the optical thickness of the medium; ii) the wall emissivities, that rule the wall reflection effects;

iii) the temperature level, that controls the non linear radiative fluxes; iv) the Reynolds number that controls the weight of turbulence transport and the balance of gas-gas and gas-wall radiative effects.

• Radiation effects on the fluctuation field within a turbulent channel flow have also been analyzed. Firstly, radiation reduces the enthalpy fluctuation and influences the turbulent transport of enthalpy fluctuation by decreasing production terms and introducing additional radiative dissipative terms. Various conditions lead to distinctive profiles of wall-scaled quantities. A new radiation-based scaling has been proposed, it improves the agreement between profiles of turbulent heat flux and other quantities in channel flows with/without accounting for radiative transfer. From this new scaling, a model for the turbulent Prandtl number is proposed to take into account radiation effects.

• In order to accurately predict the influence of radiation on the temperature profile in turbulent boundary layers, a new one-dimensional wall model for large-eddy simulation (LES) is proposed. The wall model describes the inner boundary layer which cannot be resolved by the LES. The radiative power source term is calculated from an analytical expression of the intensity field within the inner layer. Wall stress and conductive heat flux predicted by the wall model are fed back to the large-eddy simulation in the outer layer which is coupled to OERM to account for radiation.

First results have shown that: i) A semi-local scaling of the mixing-length is necessary to account for the variations of gas properties, such as mass density and dynamic viscosity, within the inner layer; ii) A non-constant turbulent Prandtl number gives better results; iii) In the studied cases, the mean incoming intensity at the outer boundary of the wall model has required the discretization of the corresponding solid angle into at least two angular sectors to accurately account for anisotropy of incident radiation field; iv) The point where information is transmitted from LES to the wall model has been moved away from the wall while remaining in the validity range of the wall-model. This is done in order to decrease the mismatch in the obtained wall laws due to numerical and sub-grid errors in the wall vicinity.

Then, the combination of the proposed wall model and LES/Monte-Carlo method has been compared to DNS/Monte-Carlo results on turbulent channel flows. Coupled results have also been compared to those obtained with a standard wall model which does not account for radiation effects within the inner layer. Results for high pressure cases show that the new wall-model greatly improves the accuracy of the predicted results compared to a standard wall model. The same trend is retrieved in atmospheric cases even though the magnitude of difference is smaller than the one at high pressure.

• Although all the cases considered in the previous study show an important radiation effect, the magnitude of radiative transfer in real applications depends on the flow condition. It could then be very useful to estimate the importance of radiation before performing an accurate simulation which is generally very expensive. In order to realize this, a criterion is proposed to assess the weight of wall radiative flux compared to wall conductive flux. Moreover, to predict whether the wall model accounting for radiation is required for a given flow, another criterion based on the change of nondimensional temperature due to radiation (scaled in wall units) is also developed. These two criteria are then validated with DNS/Monte-Carlo results. Both criteria are based on flow fields resulting from a simulation performed without radiation.

In order to carry out a parametric study, a k-ǫ model and a formula for turbulent Prandtl number are applied to simulate the mean velocity and temperature field of channel flows under different flow conditions. And the corresponding criterion values are computed and compared.

The results generalize the previous observations: An increased optical thickness of medium enhances radiation effects and higher values for both criteria are obtained as hydraulic diameter of channel or/and pressure increase. When wall emissivity changes, the trend of the two criteria are opposite since the wall radiative flux is mainly due to gas-wall radiation while change of wall-scaled temperature, on the other hand, is generally dominated by gas-gas radiation. The two criteria generally decrease with Reynolds number, while a small negative value of the criterion for wallscaled temperature can be achieved when gas-wall radiation overcomes gas-gas radiation in flows characterized by high Reynolds number.

Perspectives

The perspectives for this work include the following points:

• The radiation effect in a reactive combustion system should be investigated. Due to higher temperature fluctuations in these systems, the influence of TRI is known to be even more important.

• The turbulent Prandtl number model proposed in chapter 5 has actually been developed after the wall-model and it should be combined with the wall-modeled LES to further increase the accuracy of the approach. Its impact is small in the studied conditions where radiation was predominant but can be expected to be more noticeable when radiation and turbulent heat transfer are of the same order in magnitude.

• The wall modeled-LES can be applied to more complex systems such as combustion systems. Other physical effects, such as chemistry, pressure gradient or non-equilibrium boundary layers, could be included in the wall model. Finally, under other configurations, the requirement on the discretization of the solid angle in the wall-model might be different.

• The proposed wall model could be implemented in RANS which is widely used for engineering applications.

• The criteria maps built in chapter 7 should be used to determine operating conditions of furnaces, gas turbines, rocket engines, ... where the observed effects of radiation on the boundary layer structure must be accounted for in numerical simulations.

Publications related to the thesis

Several results obtained during this thesis have led to the following international publications: Two other publications associated with chapter 5 and 7 are planned. Additionally, several posters have been presented in international conferences:

• Y.F. Zhang, R. Vicquelin, O. Gicquel, and J. where the total heat flux q tot is defined as

q tot = λ ∂T ∂y -ρ h ′′ s v ′′ .
(A.17)

Eq. (A.17

Low Mach-number approximation

If the gas velocity is small compared to the sound speed (Mach number Ma < 0.2), the low Mach number approximation can be used where the density is decoupled from the pressure and hence acoustics is filtered out of the computation.

In order to derive the low Mach number approximation, the non-dimensional variables are expanded as power series of the Mach number. The variables are normalized as follow

Streamwise momentum source term to control the bulk Reynolds number

The momentum equation in the streamwise direction (Eq. 3.2 when i=1) can be recast as

∂ ∂t (ρu) + ∂ ∂x j (ρu i u) = - ∂p ∂X + ∂τ i1 ∂x i + S 1 (D.1)
where a homogeneous source term S 1 is added to compensate for viscous dissipation and enforces the flow to have a target bulk Reynolds number Re t b . In standard flows with constant flow properties and without multi-physical phenomena (chemistry, radiation, ...), non-dimensional equations can be written and this source term is directly specified by the Reynolds number. However, in more complex flows such as the studied ones, equations are kept in their dimensional form and S 1 must be determined differently. A first method consists in taking a fixed constant in time S ref [START_REF] Ghosh | Effects of radiative heat transfer on the structure of turbulent supersonic channel flow[END_REF][START_REF] Cabrit | Direct simulations for wall modeling of multicomponent reacting compressible turbulent flows[END_REF] that is either chosen arbitrarily, leading to an unknown and different Reynolds number in each different configuration, or more carefully evaluated from the friction coefficient using experimental/theoretical formulas, functions of the Reynolds number. Nonetheless, in the latter, the final Reynolds number that is obtained remains different from the intended one since formulas that are considered correspond to simple flows. The deviation from the target value increases as the studied flow is more and more complex. A second method consists in dynamically adapting the source term value after each iteration so that the Reynolds number is brought towards its target value [START_REF] Bocquet | A compressible wall model for large-eddy simulation with application to prediction of aerothermal quantities[END_REF][START_REF] Cabrit | Direct simulations for wall modeling of multicomponent reacting compressible turbulent flows[END_REF]. The following procedure has been reported in [START_REF] Cabrit | Modélisation des flux pariétaux sur les tuyères des moteurs à propergol solide[END_REF][START_REF] Bocquet | A compressible wall model for large-eddy simulation with application to prediction of aerothermal quantities[END_REF]): Here, another approach is proposed which ensures the channel flow to converge exactly and efficiently to the target bulk Reynolds number. The idea is to derive a second order ordinary differential equation with constant coefficients for Re b whose time response can then be controlled. Integration of Eq. D.1 over the whole computational domain V gives rise to d dt

S n+1 1 = S ref + ρ t b u t b - V ρ n u n dV ′ /V /
V ρu dV ′ = τ i1 • n dS + S 1 V (D.3)
where n is the surface normal vector. Note that the integration of pressure gradient term is null since the periodic boundary condition is applied in X directions, same for the convective flux of momentum.

The integrated term on the left side of Eq. D.3 can be expressed as The final formulation for the source term is therefore: Using a first order discretization, the time differential equation for S 1 (t) gives a new value S n 1 at each new iteration. This approach has been applied in the reported DNS and LES.

Discretization of main equations and treatment of source terms for RANS in channel flows

In section 7.2, a RANS approach is used to compute channel flows. In this part, the discretization of the main equations and the treatment of different source terms is detailed. If a second-order central difference scheme is used for the derivative terms and if the mean value of the source term is taken as the point value S Φ,i , Eq. (E.2) then writes [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF])

E.1 Discretization

D Φ,i+ 1 2 Φ i+1 -Φ i y i+1 -y i -D Φ,i-1 2 Φ i -Φ i-1 y i -y i-1 + S Φ,i ∆y i = 0. (E.3)
Generally, the source term can be expressed as (see section E.2) S Φ,i = S C Φ ,i + S P Φ ,i Φ i , (E.4) and then Eq. (E.3) is recast into the following form:

a i Φ i = a i+1 Φ i+1 + a i-1 Φ i-1 + b i , (E.
D Φ,i+ 1 2 = f i D Φ,i + 1 -f i D Φ,i+1
-1

with f i = ∆y + i ∆y i (E.10) while the calculation of D Φ,i-1 2 is similar. Equation (E.5) is solved with an iterative procedure: If the superscript n + 1 and n denote the present and previous iterative step respectively, Eq. (E.5) becomes (E.11) where i = 1, 2, ..., N and N is the number of grid points. In this study, the TriDiagonal-Matrix Algorithm (also called Thomas algorithm) [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF]) is applied to solve this linear algebraic equations (E.11).

a n i Φ n+1 i = a n i+1 Φ n+1 i+1 + a n i-1 Φ n+1 i-1 + b n i ,

E.2 Source terms

In the main equations of k-ǫ model, a source term appears in addition to other different derivative terms. For numerical stability, a proper treatment of these source terms is carried out of the iterative procedure. Moreover, the source terms of momentum and energy balance equations are also discussed since they respectively control the practical Reynolds number and temperature value on the central line of channel. Since the k-ǫ model is notoriously stiff [START_REF] Wilcox | Turbulent modeling for CFD[END_REF], spatial attention should be paid to corresponding source terms.

As detailed in Ref. [START_REF] Huang | An implicit navier-stokes code for turbulent flow modelling[END_REF], the source terms of k and ǫ can be linearized into positive part and negative part as This linearization ensures the positive value of k or ǫ and stabilizes the iterative method. Moreover, a threshold value 10 -12 is used as the minimum value of k in the iterative process. The final flow field results are not affected by this threshold value since k is always much larger than 10 -12 when converged.

S k,i = µ t

E.2.2 Source term for u

Uniform source terms S u and S T are added in the momentum and energy balance equations respectively to obtain a target practical Reynolds number Re t D h and a target central temperature T

  support and love in these years.

  c p Thermal capacity at constant pressure [J•kg -1 •K -1 ]. e Optical thickness [-]. f Probability density function [-]. h Enthalpy per unit mass [J•kg -1 ]. I Radiative intensity [W• m -2 • sr -1 ]. k Turbulent kinetic energy [m -2 • s -2 ]. volume [W•m -3 ]. q Heat flux [W•m -2 ]. Q Energy source term [W•m -3 ]. Re Reynolds number [-]. R Cumulative distribution function [-]. R Indicator for wall radiative heat flux [-]. r Distance [m]. S i Momentum source term [N•m -3 ] or surface [m 2 ]. s Propagation direction unit vector[-]. s Position of a photon bundle [m]. t Time [s]. T Temperature [K]. u Streamwise velocity component [m•s -1 ]. u i , u j Velocity vector [m•s -1 ]. V Volume [m 3 ]. X , Y , Z Cartesian coordinates [m].x i Coordinate vector (tensorial) [m]. y Distance to a wall [m].

IntroductionFigure 1 :Figure 2 :

 12 Figure 1: Energy flow diagram for one operating point of an industrial furnace, illustrating the four different contributions to output and the effect of wall-to-wall radiation exchange (from Ref. (Enomot et al. 1975)).

Figure 3 :

 3 Figure 3: Scheme of basic heat processes related to combustor liner.

Figure 4 :

 4 Figure 4: Scheme of grid system for wall-resolved LES (a) and for wall-modeled LES (b).

Introduction 21 Figure 5 :

 215 Figure 5: Instantaneous temperature fields without (top) and with (bottom) radiation (from Ref. (dos Santos et al. 2008)).

Figure 1

 1 Figure 1.1: A schematic description of radiation energy evolution.

  Figure 1.2: A schematic description of scattered beam.

τ ν = s 0 (

 0 κ ν + σ sν )ds, (1.4) the RTE can be written in a simple form asdI ν dτ ν + I ν = S ν (τ ν , s),(1.5)

Figure 1

 1 Figure 1.3: A schematic description of the reflection.

Figure 1

 1 Figure 1.4: A schematic description of different MC methods (the solid line means the radiative power is calculated directly while the dashed line indicates the computation is based on reciprocity principle).

  MONTE-CARLO METHODsum of the exchange powers P exch ij between i and all the other cells j, i.e:

  Figure 1.5: Radiative exchange between two differential volume elements.

Figure 1 . 6 :

 16 Figure 1.6: Elementary cells for reciprocal Monte-Carlo method.

  Figure 1.7: Principles of computation of ERM, ARM and FM: Cell i is the cell where radiative power is computed; The continuous arrows start from emission points and reach absorption points; The dashed arrows deals with the reciprocity principle.

  indicates that the computation of P i depends on the emission of all the cells in the enclosure (as shown in Fig.1.7). To reduce the variance of cell V i , "thesis_yufang" -2013/11/27 -19:31 -page 43 -#43

FMFigure 2 . 1 :

 21 Figure 2.1: Principles of computation of ERM, ARM and FM: Cell i is the cell where radiative power is computed; The continuous arrows start from emission points and reach absorption points; The dashed arrows deal with the reciprocity principle.

  Figure 2.2: Mean radiative power and standard deviation over N x N z points ; : reference solution; -▽-: ERM; : OERM.

Figure 2 . 3 :

 23 Figure 2.3: Spectral emitted and absorbed power for a cell of temperature 700K;:P e ν ; :P a ν .

Figure 2 . 4 :

 24 Figure 2.4: Spectral emitted and absorbed power for a cell of temperature 2500K; symbols as figure 2.3.

Figure 2 . 5 :Figure 2 . 7 :

 2527 Figure 2.5: Mean radiative power and standard deviation over N x N z points (case 1); symbols as figure 2.2.

Figure 2 . 8 :Figure 2 . 9 :

 2829 Figure 2.8: Mean radiative power and standard deviation over N x N z points (case 3); symbols as figure 2.2.

Figure 2 . 10 :

 210 Figure 2.10: Mean radiative power value and standard deviation over N x N z points (case 4); symbols as figure 2.2.

Figure 3 . 1 :

 31 Figure 3.1: Typical energy spectrum of turbulence (k c is the cut-off wave number in LES; k d = 2π/η, where η is the Kolmogorov length scale).

  3.When dealing with numerical simulations of turbulent flows, three different approaches are generally considered: Direct Numerical Simulation, Large Eddy Simulations and Reynolds Averaged Navier-Stokes Equations. These three approaches are discussed in the following. Models and numerical approaches for the flow 63

  Models and numerical approaches for the flow 65

  Models and numerical approaches for the flow 67

  Models and numerical approaches for the flow 71 with A + = 25.

  introduced stochastic forcing in the interface region of RANS and LES which could accelerate the "thesis_yufang" -2013/11/27 -19:31 -page 72 -#72

  Models and numerical approaches for the flow 73

  .36) By analyzing the importance of the two terms on the left side of Eq. (3.36), the turbulent boundary layer can be divided into two main parts: a viscous sublayer and a logarithmic sub-layer[START_REF] Schlichting | Boundary layer theory[END_REF][START_REF] Pope | Turbulent Flows[END_REF] (as shown in Fig.3.2).

  Figure 3.2: Velocity profile of turbulent boundary layer in wall units.

  .56) "thesis_yufang" -2013/11/27 -19:31 -page 77 -#77 Models and numerical approaches for the flow 77 Equation (3.34) can be transformed into du * dy * + (κy * ) 2 du * dy *

"

  thesis_yufang" -2013/11/27 -19:31 -page 78 -#78

  .68) "thesis_yufang" -2013/11/27 -19:31 -page 79 -#79

  YALES2 employs several iterative solvers to solve the Poisson equation. They are: Preconditioned Conjugated Gradient (PCG), Bi-Conjugate Gradient stabilized (BICGSTAB), Bi-Conjugate Gradient stabilized(2) (BICGSTAB2), Deflated PCG (DPCG) and Deflated BiCGstab(2) (DBICGSTAB2). Among "thesis_yufang" -2013/11/27 -19:31 -page 80 -#80 COUPLING BETWEEN YALES2 AND THE MONTE-CARLO METHOD them, DPCG is used in this work with a convergence criterion of 1E-8. The mass conservation equation is implicitly fulfilled by solving the Poisson equation.

"Figure 3 . 3 :

 33 Figure 3.3: Schema of coupling.

Figure 3

 3 Figure 3.4: Instantaneous temperature profile along a line at different number of iterations (position of the line x = 0.0, z = 0.0).

  "thesis_yufang" -2013/11/27 -19:31 -page 83 -#83

  DNS of turbulent channel flow: effect of radiation on the mean flow fields 85

  ; dos[START_REF] Dos Santos | Coupled large eddy simulations of turbulent combustion and radiative heat transfer[END_REF] have shown an impact of radiation on turbulent flame dynamics;[START_REF] Ghosh | Effects of radiative heat transfer on the structure of turbulent supersonic channel flow[END_REF] have studied the effects of radiation in a turbulent compressible channel flow and have shown that radiation modified the Reynolds stresses in the near wall layer even for an optically thin medium. Using a low-Reynolds RANS model,[START_REF] Soufiani | Radiation effects on turbulent heat transfer in channel flows of infrared active gases[END_REF] have also carried out a channel flow simulation where the wall conductive heat flux has been found to be significantly affected by radiation while, in a similar case with DNS here,[START_REF] Amaya | Unsteady coupling of Navier-Stokes and radiative heat transfer solvers applied to an anisothermal multicomponent turbulent channel flow[END_REF] have found the effect of radiation in the boundary layer to be weak. Using a grey gas model and varying the medium optical thickness,[START_REF] Gupta | Large-Eddy Simulation of Turbulence-Radiation Interactions in a Turbulent Planar Channel Flow[END_REF] have shown that the temperature profile in a turbulent channel flow can be modified by radiation."thesis_yufang" -2013/11/27 -19:31 -page 86 -#86

  .1) "thesis_yufang" -2013/11/27 -19:31 -page 87 -#87 DNS of turbulent channel flow: effect of radiation on the mean flow fields 87

Figure 4 . 1 :

 41 Figure 4.1: Computational domain of channel flow cases. X, Y and Z are the streamwise, wall normal and spanwise directions. L X , L Y and L Z are the dimensions of the channel case in each direction. δ is the channel half-width. The lower wall (resp. upper wall) is at temperature T w,c (resp. T w,h ; T w,h ≥ T w,c ). Periodic boundary conditions are applied along X and Z.

  DNS of turbulent channel flow: effect of radiation on the mean flow fields 89 Re bReD h Tw,c(K) T w,h (K) p(atm)

Figure 4 . 2 :

 42 Figure 4.2: Mean profiles of non-dimensional streamwise velocity (a): (Kim et al. 1987); and non-dimensional temperature in wall units (b): △(Pr = 0.1), ▽(Pr = 0.71), (Pr = 2) (Kim and Moin 1987); : present results.

Figure 4 . 3 :

 43 Figure 4.3: Comparison of mean radiative power of C1R1, C1R1_GG and C1R1_GW ( : C1R1; : C1R1_GW; : C1R1_GG; color : Cold side; color : Hot side; error bars represent the standard deviation).

Figure 4 . 4 :

 44 Figure 4.4: Mean heat flux distribution for C2 (a), C2R3 (b) and comparison of these two cases (c) ( : Turbulent convective heat flux; : Conductive heat flux; • • • : Radiative energy flux; • : total heat flux; black lines: C2; gray lines: C2R3).

Figure 4 . 5 :

 45 Figure 4.5: Mean temperature profiles in global coordinates (a) and in wall units (b): ▽ DNS results (Pr= 0.71) (Kim and Moin 1987) ; : C1; : C1R1_GW; : C1R1_GG; : C1R1; in (b), : Cold side; color : Hot side.

  DNS of turbulent channel flow: effect of radiation on the mean flow fields 97

Figure 4 . 6 :

 46 Figure 4.6: Comparison of mean temperature profile (a), T + profile (b,c) of C2, C2R1_GW, C2R1_GG and C2R1 (▽: DNS results of a passive scalar of Kim and Moin : (Kim et al. 1987) (with Pr=0.71); : C2; : C2R1_GW; : C2R1_GG; : C2R1; in (b), color : Cold side; color : Hot side).

Figure 4 . 7 :

 47 Figure 4.7: Mean profile of temperature (a), T + on the cold side (b) for C2, C2R1, C2R2 and C2R3 ( ▽: DNS results of a passive scalar of Kim and Moin (Kim et al. 1987) (with Pr=0.71); : C2; : C2R1; • : C2R2; • • • : C2R3).

  4.3 ) corresponds to two extreme cases: Very reflecting cold wall of emissivity 0.1 and black hot wall. The previous effects on the averaged temperature field are here amplified, as shown in Fig. "thesis_yufang" -2013/11/27 -19:31 -page 100 -#100

Figure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Mean profile of total radiative power (a), gas-gas radiative power (b), gas-wall radiative power (c) and gas-cold wall and gas-hot wall radiative power (d) of C2R1, C2R2 and C2R3 ( : C2R1; • : C2R2; • • • : C2R3; in (a) and (b), color : Cold side; color : Hot side; error bars in (a) represents standard deviation).

Figure 4

 4 Figure 4.10: Mean profile of radiative power for the channel hot side for C2R3 (a) and C2R4 (b) ( ••: Total; * : Gas-wall; ⋄ : Gas-gas).

Figure 4 . 11 :

 411 Figure 4.11: Mean profile of velocity (a), u + (b) and u + V D (c) of C3, C3R1 and C1 ( • : Wall function of Ref. (Kim et al. 1987) ; : C3; •• : C3R1; : C1; in (b) and (c), color : Cold side; color : Hot side).

Figure 4 .

 4 Figure 4.12: Mean profile of mean temperature(a) and T + (b) of C3, C3R1, C1 and C1R1 ( • : Wall function of Ref. (Kasagi et al. 1992) ; : C3; •• : C3R1; : C1; : C1R1; in (b), color : Cold side; color : Hot side ).

  Fig. 4.15 where they are compared to results "thesis_yufang" -2013/11/27 -19:31 -page 105 -#105

Figure 4 .

 4 Figure 4.13: Mean temperature profile (a) and T + of C4, C4R1, C1 and C1R1 (b) ( : C4; •• : C4R1; : C1; : C1R1; in (b), color : Cold side; color : Hot side).

Figure 4 . 14 :

 414 Figure 4.14: Comparison of turbulent convective heat flux (a) and radiative energy flux (b) of C4R1 and C1R1 ( •• : C4R1; : C1R1).

Figure 4 . 15 :

 415 Figure 4.15: Mean radiative power profile of C4R1 and C1R1 on the cold (left) and hot (right) sides ( : Total; : Gas-wall; : Gas-gas; black lines : C4R1; gray lines : C1R1; the error bars represent the standard deviation).

  channel flow: effect of radiation on second order moments This chapter focuses on radiation effects on fluctuation field in turbulent channel flow. The influence of radiation on enthalpy root-meansquare, on turbulent heat flux and on the budget of the corresponding transport equation are analyzed. To improve the agreement of results of non-dimensional turbulent quantities between flows without and with radiation, a new radiation-based scaling is proposed. The influence of radiation on turbulent Prandtl number is also presented and a model based on the new scaling is proposed. The details of this study are presented in the format of a paper entitled: Analysis of direct numerical simulations coupled to radiation in turbulent channel flow.

Figure 5 . 1 :

 51 Figure 5.1: Computational domain of channel flow cases with the half-width δ = 0.1 m. The lower wall (resp. upper wall) is at temperature T w,c (resp. T w,h ; T w,h ≥ T w,c ). Re b Tw,c [K] T w,h [K] p [atm] A 5850 950 1150 40.0 B 950 2050 40.0 C 11750 950 1150 40.0

Figure 5 . 2 :

 52 Figure 5.2: Profiles of mean temperature scaled by wall temperatures (a) and in wall units (b) on the cold side (blue color) and hot side (red color) in cases A (thin line) and A_R (thick line).

Figure 5 . 3 :

 53 Figure 5.3: Profiles of enthalpy root-mean-square scaled by center temperature and thermal capacity (a) and in wall units (b) on the cold side (blue color) and hot side (red color) in cases A (thin line) and A_R (thick line).

Figure 5 . 4 :

 54 Figure 5.4: Budget of enthalpy variance (cold side only) in cases A (a) and A_R (b): Production (plain line); Molecular dissipation (dashed line); Radiative dissipation (dashed-dashed-dotted line); Turbulent diffusion (dashed-dotted line); Molecular diffusion (dotted line); Density-enthalpy correlation term (dashed-dotted-dotted line).

Figure 5 . 5 :

 55 Figure 5.5: Profiles of wall-normal turbulent heat flux scaled by bulk variables (a) and in wall-units (b) on the cold side (blue color) and hot side (red color) in cases A (thin line) and A_R (thick line).

Figure 5 . 6 :

 56 Figure 5.6: Profiles of correlation between enthalpy and wall-normal velocity fluctuations in cases A (thin line) and A_R (thick line).

Figure 5 . 7 :

 57 Figure 5.7: Budget of wall-normal turbulent heat flux (cold side only) in cases A (a) and A_R (b): Production (plain line); Molecular diffusion (dotted line); Turbulent diffusion (dashed-dotted line); Pressure term (dashed-dotted-dotted line); Molecular dissipation (dashed line); Velocity-Radiative power correlation term (dashed-dasheddotted line).

Figure 5 . 8 :

 58 Figure 5.8: Profiles of the P R (a) and P R rms (b) in case A_R on the cold side (blue color) and hot side (red color). Circles : mean radiative power computed from the mean temperature field.

Figure 5 . 9 :

 59 Figure5.9: Correlation coefficients between radiative power and enthalpy (a) and between radiative power and wall-normal velocity (b) in case A_R on the cold side (blue color) and hot side (red color).

C

  Figure 5.10: Correlation coefficient defined in Eq. 5.25 in case A_R.

Figure 5

 5 Figure 5.11: Scheme of coupling effects in turbulent channel flow with radiation. The fields of velocity (u, v, w), temperature (T ) and radiation are split into mean and fluctuating components. Several quantities are considered for the radiation field: The radiative power (P R ), the spectral intensity (I ν ) and absorption coefficient (κ ν ). Main effects are represented by thick plain arrows while thin dotted arrows indicate negligible/null effects.

Figure 5

 5 Figure 5.12: Profiles of mean temperature scaled by wall temperatures (a) and in wall units (b) on the cold side (blue color) and hot side (red color) in cases A (thin plain line), A_R (thick plain line), B (thin dashed line) and B_R (thick dashed line).

Figure 5

 5 Figure 5.13: Profiles of turbulent kinetic energy on the cold side (blue color) and hot side (red color) in cases A (thin plain line), A_R (thick plain line), B (thin dashed line) and B_R (thick dashed line). (a) With standard wall-scaling. (b) With semi-local scaling.

Figure 5 . 14 :

 514 Figure 5.14: Profiles of turbulent heat flux (a) and production (b), pressure term (c) and correlation between velocity and radiative power (d) in semi-local scaled units on the cold side (blue color) and hot side (red color) in cases A (thin plain line), A_R (thick plain line), B (thin dashed line) and B_R (thick dashed line).

  from the one using standard wall-scaling. Profiles of ρv ′′ h ′′ * are compared inFig. 5.14 (a). Good agreement between results of A and B is obtained. However, cases with radiation A_R and B_R are quite different nor do they match each other. Radiation effects are strongly dependent on tem-"thesis_yufang" -2013/11/27 -19:31 -page 128 -#128

Figure 5 . 15 :

 515 Figure 5.15: Profiles of turbulent heat flux (a) and production (b), pressure term (c) and correlation between velocity and radiative power (d) using radiation-based scaling on the cold side (blue color) and hot side (red color) in cases A (thin plain line), A_R (thick plain line), B (thin dashed line) and B_R (thick dashed line).

Figure 5 . 16 :

 516 Figure 5.16: Profiles of enthalpy rms (a) and production (b), molecular dissipation (c) and radiative dissipation (d) using semi-local scaling on the cold side (blue color) and hot side (red color) in cases A (thin plain line), A_R (thick plain line), B (thin dashed line) and B_R (thick dashed line).

Figure 5 . 17 :

 517 Figure 5.17: Profiles of enthalpy rms (a) and production (b), molecular dissipation (c) and radiative dissipation (d) using radiation-based scaling on the cold side (blue color) and hot side (red color) in cases A (thin plain line), A_R (thick plain line), B (thin dashed line) and B_R (thick dashed line).

Figure 5 .Figure 5

 55 Figure 5.19 compares cases A and A_R to cases C and C_R, characterized by a larger Reynolds number, interns of mean temperature profiles. On the one hand, radiation effects are reduced by the increased influence of turbulent

Figure 5 . 19 :Figure 5 . 20 :Figure 5 . 21 :Figure 5 . 22 :Figure 5 . 23 :

 519520521522523 Figure 5.19: Profiles of mean temperature scaled by wall temperatures (a) and in wall units (b) on the cold side (blue color) and hot side (red color) in cases A (thin plain line), A_R (thick plain line), C (thin dashed-dotted line) and C_R (thick dasheddotted line).

  Figure 5.24: Profiles of the turbulent Prandtl number on the cold side (blue color) and hot side (red color) in cases A (thin plain line), A_R (thick plain line), B (thin dashed line), B_R (thick dashed line), C (thin dashed-dotted line) and C_R (thick dashed-dotted line).

Figure 5

 5 Figure 5.25: Profiles of time scale ratio R in cases A (thin plain line) and A_R. In case A_R, different definitions of the thermal time scale are based on: Molecular dissipation only (thick plain line); Radiative dissipation only (thick dashed line); Both dissipation terms (thick dashed-dotted line).

Figure 5

 5 Figure 5.27: Profiles of modeled turbulent Prandtl number on the cold side (a) and hot side (b) compared to DNS in cases A_R (DNS: Circles; Model: Plain line) , B_R (DNS: Down-pointing triangles; Model: Dashed line) and C_R (DNS: Up-pointing triangles; Model: Dotted line).

  DNS of turbulent channel flow: effect of radiation on second order moments 145

  Figure 5.28: Relative variation of temperature T max rms /T (y + BL ) as a function of the Reynolds number for several values of T b /T w .

  for LES accounting for radiation effectsAs detailed in chapter 4, temperature profile in turbulent boundary layer and the corresponding wall conductive heat flux are strongly modified by radiation. In order to take into account this effect of radiation, a new wall model for large eddy simulation (LES) is proposed in this chapter. After a separate validation of different components of the proposed wall model, full-coupled wall-model LES/Monte-Carlo simulation of turbulent channel flow is performed. And the obtained results are compared with DNS/Monte-Carlo results from chapter 4 to assess the accuracy of the new model. This study is detailed in the following submitted publication: Y.F. Zhang, R. Vicquelin, O. Gicquel and J. Taine (submitted). A wall model for LES accounting for radiation effects. International Journal of Heat and mass Transfer (2013).

  in the turbulent transport model. In other wall models "thesis_yufang" -2013/11/27 -19:31 -page 151 -#151 A wall model for LES accounting for radiation effects 151

  are considered to validate the proposed "thesis_yufang" -2013/11/27 -19:31 -page 152 -#152 WALL-MODELED LES COUPLED TO RADIATION

0 Figure 6 . 1 :

 061 Figure 6.1: Scheme of wall model grids embedded within the LES grid. In this example, the wall distance y w for the coupling point corresponds to the first LES grid point off the wall.

  A wall model for LES accounting for radiation effects 155

Figure 6 . 2 :

 62 Figure 6.2: Definition of forward and backward intensities I + ν and I - ν (µ 1 and µ 2 are the cosine of polar angle θ 1 and θ 2 respectively; y = 0 and y = y w are the two boundaries of the 1D model).

  WALL-MODELED LES COUPLED TO RADIATION is assumed isotropic. At any point y of the inner layer, the radiative power per unit volume P R (y), difference between absorbed and emitted powers, writes P R (y) = P a (y) -P e (y)

τ

  νy ′ y (µ) = exp[-e ν (y ′ , y)/µ] with : e ν (y ′ , y) = y y ′ κ ν (y ′′ )dy ′′ . (6.19)

0 κ 0 κ

 00 A wall model for LES accounting for radiation effects 157 and the generalized incomplete function defined by E n (X, [µ j , µ j+1 ]) = of the radiative power P R (y) is simply given byP R (y) ≈ 2π e ν (y, y w ), [µ j , µ j+1 ]) I ν (y w , [µ j , µ j+1 ]) ν (y ′ )I • ν (y ′ ) E 1 (e ν (y ′ , y)) dy ′ dν, ν (y) I • ν (y) dν. (6.22)where N µ is the number of angular sectors [µ j , µ j+1 ] used to discretize the 2π steradians associated with the incoming intensity at a point y w . I ν (y w , [µ j , µ j+1 ]) is the average value of I - ν (y w , µ) over the range [µ j , µ j+1 ] 1 . The mean anisotropic spectral incoming intensity field I ν (y w , [µ j , µ j+1 ]) at any grid point y w is determined by the Monte Carlo method, as detailed in Appendix.

  SEPARATE VALIDATION OF THE TURBULENCE AND RADIATION WALL MODEL COMPONENTScouple LES, Monte-Carlo model and the radiation and turbulence wall models is schematized in Fig.6.3.

Figure 6 . 5 :

 65 Figure 6.5: Mean velocity profile in wall units, i.e. {u}/u τ , of case B (only the cold side) and Ref. (Hoyas and Jimenez 2006) (circles: DNS reference data; black line: Semi-local mixing-length model; gray line: Standard wall mixing-length model).For the sake of clarity, curves of the different cases are translated.

Figure 6 . 6 :Figure 6 . 7 :Figure 6 . 8 :

 666768 Figure 6.6: Mean temperature profile in wall units, i.e. |{T } -T w |/T τ with T τ = |q cd w |/( ρ w c pw u τ ), on the cold side of case B (circles: DNS datas from Ref. (Zhang et al. 2013); black lines: Semi-local mixing-length model; gray lines: Standard mixinglength model; dashed-dotted line: Constant Pr t ; dashed line: Eq. (6.5) for Pr t ; plain line: Eq. (6.6) for Pr t .

  A wall model for LES accounting for radiation effects 163 associated with N µ = 1 are practically acceptable.

Figure 6 . 9 :Figure 6 . 10 :

 69610 Figure 6.9: Radiative power on the cold side of case B with y w = 0.05δ (a), y w = 0.1δ (b) and y w = 0.2δ (c) (circles: DNS datas from Ref. (Zhang et al. 2013); dashed-dotted line: N µ = 1; dashed line: N µ = 2; plain line: N µ = 3).

Figure 6 . 11 :Figure 6 . 12 :Figure 6 . 13 :

 611612613 Figure 6.11: Rms of velocity components u (a), v (b), w (c) and temperature (d) in wall units on the cold side of case A (circles: DNS datas from (Zhang et al. 2013); dashed-dotted line: Y1; dashed line: Y2; plain line: Y3).

Figure 6 . 14 :

 614 Figure 6.14: Mean velocity profile in wall units on the cold side (a), mean temperature profile in wall units on the cold side (b) and hot side (c) of cases D_RAD03 and D_RAD08 (circles: DNS datas from (Zhang et al. 2013); dashed-dotted line: No wall model; dashed line: Standard wall model; plain line: New wall model). Curves of the different cases are translated.

Figure 6

 6 Figure 6.15: Scheme showing how the mean incoming intensityI ν (y w , [µ j , µ j+1 ]),boundary condition of the wall model, is computed from the reciprocal Monte-Carlo method by using a virtual wall. µ j and µ j+1 correspond to the cosine of angles θ j and θ j+1 respectively.

Figure 6 . 16 :

 616 Figure 6.16: Profiles of turbulent Prandtl number (a,c) and scaled mean temperature profile (b,d) on the cold (a,b) and hot (c,d) sides in case B_R (circles: DNS; black dashed line: New turbulent Prandtl number model Eq. (5.45); gray plain line: Eq. (6.6) for turbulent Prandtl number).

  Practical criteria to compare radiative and conductive fluxes and their coupling effectsA criterion controlling the weight of radiative flux on walls is proposed to assess the necessity of performing radiation simulation under a given flow condition. Similarly, to predict whether the wall model accounting for radiation is required for a given flow, another criterion based on the change of non-dimensional temperature due to radiation (scaled in wall unit) is also developed. The validity of these two criteria is confirmed by an analysis of fully-coupled DNS/Monte-Carlo results. Mean flow fields of many turbulent channel flows at different flow conditions are then resolved with a k-ǫ model and a formula for turbulent Prandtl number.

  Having a good estimation of wall heat flux is mandatory in engineering applications. This flux can be decomposed into two parts: The conductive heat flux and the radiative flux. To evaluate the importance of wall radiative flux com-

  .11) Here the term on the left hand side, varying with position y, is the relative change of wall-scaled temperature due to radiation. The position y = "thesis_yufang" -2013/11/27 -19:31 -page 178 -#178 is here chosen to quantify the importance of this change and a criterion ∆ ref t + ,coupled y=0.1δ

Figure 7 . 1 :

 71 Figure 7.1: Results of criteria for wall radiative flux (a) and for T + (b).

Figure 7 . 2 :

 72 Figure 7.2: Velocity profiles (a) and temperature profiles (b) of the case of Kim and Moin(Kim et al. 1987).

Fig

  Fig. 7.3 (a) shows the contour plot of criterion R q when practical Reynolds number Re D h and hydraulic diameter of channel vary. The wall emissivity is 0.8 and the pressure is set to 40 atm in order to have relatively optically thick medium. a criterion value of 0.2 is chosen as a threshold to start a radiation calculation since a weight of 20% of radiative flux is already important. As expected in such a case, the region where R q is higher than 0.2 and hence where a radiation simulation is needed, covers the main part of the domain (gray zone in the figure). Moreover, for a certain Reynolds number, the importance of wall radiative flux increases with hydraulic diameter of channel due to an increased optical thickness. On the other hand, weight of wall radiative flux becomes

Figure 7

 7 Fig. 7.3 (a) shows the contour plot of criterion R q when practical Reynolds number Re D h and hydraulic diameter of channel vary. The wall emissivity is 0.8 and the pressure is set to 40 atm in order to have relatively optically thick medium. a criterion value of 0.2 is chosen as a threshold to start a radiation calculation since a weight of 20% of radiative flux is already important. As expected in such a case, the region where R q is higher than 0.2 and hence where a radiation simulation is needed, covers the main part of the domain (gray zone in the figure). Moreover, for a certain Reynolds number, the importance of wall radiative flux increases with hydraulic diameter of channel due to an increased optical thickness. On the other hand, weight of wall radiative flux becomes

Figure 7 Figure 7

 77 Figure 7.4: contour plots of criteria R q for heat fluxes (a) and ∆ ref t + y=0.1δ for T + (b) when change practical Reynolds number (Re D h ) and wall emissivity (pressure=40 atm; hydraulic diameter=0.4 m).

Figure 7

 7 Figure 7.6: contour plots of criteria R q for heat fluxes (a) and ∆ ref t + y=0.1δ for T + (b) when change practical Reynolds number (Re D h ) and hydraulic diameter of channel (pressure=1 atm; wall emissivity=0.8).

Figure 7 Figure 7

 77 Figure 7.7: contour plots of criteria R q for heat fluxes (a) and ∆ ref t + y=0.1δ for T + (b) when change practical Reynolds number (Re D h ) and wall emissivity (pressure=1 atm; hydraulic diameter=0.4 m).

Figure 7 Figure 7

 77 Figure 7.9: contour plots of criteria R q for heat fluxes (a) and ∆ ref t + y=0.1δ for T + (b) when change practical Reynolds number (Re D h ) and pressure (hydraulic diameter=0.4 m; wall emissivity=0.8).

Figure 7

 7 Figure 7.11: contour plots of criteria R q for heat fluxes (a) and ∆ ref t + y=0.1δ for T + (b) when change hydraulic diameter of channel and pressure (wall emissivity=0.8; practical Reynolds number =48000).

•

  Y.F. Zhang, O. Gicquel and J. Taine (2012). Optimized Emission-based Reciprocity Monte Carlo Method to speed up computation in complex systems. International Journal of Heat and mass Transfer 55 (2012) 8172-8177; • Y.F. Zhang, R. Vicquelin, O. Gicquel and J. Taine (2012). Physical study of radiation effects on the boundary layer structure in a turbulent channel flow. International Journal of Heat and mass Transfer 61 (2013F. Zhang, R. Vicquelin, O. Gicquel and J. Taine. A wall model for LES accounting for radiation effects. submitted to International Journal of Heat and mass Transfer (2013);

and v = 0 .

 0 h ′′ ) ∂z . (A.13)Here, for periodic channel flows, one has ∂(ρ u ′′ h ′′ ) Hence the terms on the left side is then reduced to only one term:∂(ρ v ′′ h ′′ )∂y . Furthermore, the first three terms on the right hand can be transformed as

  ) is the equilibrium energy equation for turbulent boundary layers in channel

x = x 0 x 0 +

 00 * , v = v 0 v * , t = x 0 /v 0 t * , ρ = ρ 0 ρ * , T = T 0 T * , h = v 2 0 h * , p = ρ 0 r 0 T 0 p * (B.1)where the subscript 0 denotes a reference value of the corresponding quantity. Substituting these equations into the governing equations yields[START_REF] Julien | Simulation numérique d'écoulements réactifs instationnaires à faibles nombres de Mach[END_REF] γ is the heat capacity ratio and Re is the Reynolds number. Considering that the Mach number is small, the non-dimensional pressure can be expressed "thesis_yufang" -2013/11/27 -19:31 -page 200 -#200 Ma p * 1 + γMa 2 p * 2 + O(Ma 3 ), (B.6) The set of equations for low-Mach number flows is obtained by i) Substituting this equation into Eqs. (B.3) and (B.4); ii) Identifying systems of equations for each order in the Mach number; iii) and keeping only the dominant orders for low Mach-number, i.e. Ma -2 , Ma -1 and Ma 0 . The obtained system of equations B.7) indicates that the pressure p * 0 = p * 0 (t) is uniform in space while Eq. (B.10) implies that p * 1 has the same feature as p * 0 . Therefore, it is reasonable to set p * 1 = 0 and hence from Eq. (B.6), the pressure now can be decomposed into two parts: a constant pressure p * 0 and the pressure fluctuation γMa 2 p

  Figure C.1: Two-point correlations in streamwise and spanwise direction for C2 at different wall normal distance.

V 0

 0 ρu dV ′ = ρ b u b V = µ b Re b V δ (D.4)where δ is the half-height of the channel and subscript b is related to the bulk quantities as aforementioned. Moreover, we haveτ i1 • n ds = -(|τ w,c | + |τ w,h |)S w (D.5)where τ w,c and τ w,h is the friction stress at the cold and hot wall respectively and S w is the surface area of each wall of the computational domain. For simplicity, the change of gas properties are not considered in the calculation of wall friction stress, which yields:|τ w,c | = |τ w,h | = |τ w |. Meanwhile,the mean wall friction stress |τ w | = (|τ w,c | + |τ w,h |)/2 can be estimated from the friction coefficient C f as |τ w | = ρ is approximately evaluated from the practical Reynolds number Re b,D h as (W.M. Kays and Weigand 2004) C f = 0.046 Re -0.2 D h for 10 4 < Re D h < 5 × 10 6 (DThe solution of Eq. D.15 writes Re b = Re t b + e αt [C 1 cos(βt) + C 2 sin(βt)] (D.16)with α = -p/2, β = 4q -p 2 /2. This damped regime is not desirable since the computational time to reach the target value is penalized by oscillations.Critically damped and overdamped regime:p 2 -4q ≥ 0The solution of Eq. D.15 writesRe b = Re t b + C 1 e r 1 t + C 2 e r 2 t , for p 2 -4q > 0 Re b = Re t b + (C 1 + C 2 t ) e r 1 t , for p 2 -4q = 0 (D.17) with r 1 = -p + p 2 -4q 2 , r 2 = -p -p 2 -4q2. (D.18) Among both regimes, the critical regime (p 2 -4q = 0) gives the shortest delay to reach the target value (see Fig. D.1). The corresponding optimized value of a 1 in Eq. D.11 is then

  Figure D.1: Change of Re b with number of iterations when ratio of p 2 /4q varies (the target Reynolds number is 6000; plain line: p 2 /4q = 1; dashed line: 10; dashed-dasheddotted line: 50; dashed-dotted line: 100; dotted line: 500).

Figure

  Figure E.1: A schematic of discretization.

  Discretization of main equations and treatment of source terms for RANS in channel flows 211 E.2.1 Source terms for k and ǫ

tc.

  When the flow is stabilized at Re D h , the integral value of S u over the whole channel flow should equal the sum of shear stresses on the two walls. Since the velocity profile is symmetric in this study, the shear stress on the two wall are the same. Therefore,S u = τ w /δ, (E.16)where the wall shear stress τ w can be calculated from the skin friction coefficient C , indice b relates to bulk quantities and C f can be approximated as (W.M.[START_REF] Kays | Convective Heat and Mass Transfer[END_REF])C f = 0.046Re -0.2 D hfor 3 × 10 4 < Re D h < 10 6 (Eis the hydraulic diameter and it equals to 4δ for a channel flow. Therefore, the source term for velocity expressesD h ) 1.8 -(Re n D h ) 1.8 . (E.21)E.2.3 Source term for TIf the temperature at two walls are equal, a source term is also needed for compensating heat losses at walls and for sustaining the central temperature value T t c . This source term S T then writesS T = -q cd w /δ. (E.22)where wall conductive heat flux q cd w can be calculated from the Nusselt number Nu, which is defined temperature T b and hydraulic diameter D h are chosen as the reference temperature and length respectively. Similarly, the formula from Ref.(W.M.[START_REF] Kays | Convective Heat and Mass Transfer[END_REF]) is applied for Nu:Nu = 0.023Re 0.8 D h Pr 0.88 + 2.03 Pr 2/3 -0.78 Re -0.1 D h (E.24)with the valid rang of 3 × 10 4 < Re D h < 10 6 and 0.6 < Pr < 6. Then the source term writes S T = Nu 4λ b (T b -T w ) law profile for mean temperature: T b ≈ 7 8 (T c -T w ) + T w , the source term S T at step n + 1 can then be calculated as:S n+1T,i = S n T,i + Nu 32λ b

Table 2

 2

			.1: Validation cases
			case 1 case 2 case 3 case 4
	L(m)		0.2	0.2	4	4
	ǫw		1.0	1.0	1.0	1.0
	Tc(K)		2500	500	2500	500
	Tw(K)	500	2500	500	2500
	Table 2.2: Averaged wall flux over two walls, and in parentheses, the corresponding
	standard deviation, the unit of wall flux is kW/m 2
			case 1	case 2	case 3	case 4
	reference solution		54.00	-100.5	195.6	-317.2
	ERM	52.69( 33.1)	-99.7(23.5)	174.9(185)	-316.0(36.4)
	OERM		53.39(12.2)	-101.8(23.4)	193.7(22.0)	-315.4(38.7)
	FM(Tessé et al. 2002)		53.98(3.02)	-		194.7(17.7)	-
	ERM(Tessé et al. 2002) 52.46(39.93)	-	182.5(283.7)	-
	ARM(Tessé et al. 2002)		53.98(2.98)	-		194.8(17.6)	-
	2.5 Results for real gases and discussion
	OERM results are characterized in this section from four typical cases which
	involve a CO 2 -H 2 O-N 2 mixture, at 1 atm, of molar fractions: x CO 2 = 0.116,
	x H 2 O = 0.155. The same system and temperature field as in Sec.2.4.1 are
	considered. The other conditions are shown in table 1. Gas spectral radiative
	properties are treated by a CK model

Table 4 .

 4 1: Channel flow parameters. In case C0, non-dimensional passive scalars are considered. Reynolds numbers are given with an accuracy of 0.5 %.

  Table 4.2: Discretization of the simulated cases: n X , respectively n Y and n Z , is the number of points in the X direction, respectively Y and Z direction. Domain size: (4πδ, 2δ, 2πδ) for C0 and (2πδ, 2δ, πδ) for C1-C4. ∆X + and ∆Z + are given at the cold and hot sides for C1-C4.

	✐					✐
		"thesis_yufang" -2013/11/27 -19:31 -page 90 -#90
	✐					✐
	90	4.3. MODELS AND NUMERICAL APPROACHES
			nX × nY × nZ	∆X +	∆Y +	∆Z +
				hot/cold		hot/cold
		C0	120 × 100 × 120 17.85/17.85 [0.8-8.0] 8.92/8.92
		C1,C2 110 × 135 × 110	16.9/21.3	[0.8-8.0]	8.4/10.6
		C3	160 × 163 × 160	8.9/22.4	[0.8-8.0]	4.4/11.2
		C4	200 × 230 × 200	17.3/21.8	[0.8-8.0]	8.6/10.8
	✐					✐
	✐					✐

Table 4

 4 For all cases considered in the paper, a non-reacting CO 2 -H 2 O-N 2 gas mixture flows, in developed turbulent regime, through a plane channel in different conditions. The molar fractions of CO 2 , H 2 O and N 2 are 0.116, 0.155 and 0.729. Cases involving radiation cannot be entirely characterized by non-dimensional numbers, contrary to cases without radiation. The value of δ, the channel half-width, is then given here: δ = 0.1 m. The dynamic viscosity µ is computed as a function of temperature from the CHEMKIN package[START_REF] Kee | A fortran computer code package for the evaluation of gas-phase, multicomponent transport properties[END_REF][START_REF] Kee | CHEMKIN-II: A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase Chemical Kinetics[END_REF]) for the chosen mixture composition. The thermal conductivity λ is computed from a Prandtl number Pr. In all considered configurations, the Prandtl number is very close to the chosen value 0.71. The flow computational cases, called Cn (n= 1 to 4), are defined in Table4.1 by a set of bulk Reynolds number Re b , pressure and wall temperatures. Similarly, cases which include radiative energy transfer are called Rm (m= 1 to 4) and are defined by the emissivities (ε 1 , ε 2 ) of the opaque walls given in Table4.3. Consequently, a computation case without radiation in conditions n is called Cn and a computation that accounts for radiation effects, in conditions n and m, CnRm. The associated spatial discretizations are defined in Table4.2. Gas radiative properties are treated in a correlated manner by the CK approach

	.3: Wall emissivities in radiative conditions R1, R2, R3 and R4.
	4.3.5 Coupled simulation

Table 5

 5

.1: Channel flow parameters: Bulk Reynolds number Re b , wall temperatures and pressure. Cases without radiation A, B and C correspond to cases C1, C3 and C4 in Ref.

[START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]

) respectively, and radiative cases A_RAD, B_RAD and C_RAD to cases C1R1, C3R1 and C4R1, respectively.

  .228( µ t / µ ) -0.0441( µ t / µ ) 2 (1 -exp(

					✐
		"thesis_yufang" -2013/11/27 -19:31 -page 154 -#154
					✐
	154	6.3. WALL-MODELED LES COUPLED TO RADIATION
	layer, two different formula proposed in Ref. (Kays 1994), both fitted from
	experimental data, are considered
	Pr t =	2.0 Pr ( µ t / µ )	+ 0.85,	(6.5)
	and
	Pr t =	0.5882 + 0	1
	✐				✐
	✐				✐

Table 6

 6 Computational domain of channel flow cases. X, Y and Z are the streamwise, wall normal and spanwise directions. L X , L Y and L Z are the dimensions of the channel case in each direction. δ is the channel half-width. The lower wall (resp. upper wall) is at temperature T w,c (resp. T w,h ; T w,h ≥ T w,c ). Periodic boundary conditions are applied along X and Z.

	.1. For cases A, B, C and D, the temperature of the two

Table 6 .

 6 1: Channel flow parameters: Bulk Reynolds number Re b , Reynolds number based on the hydraulic diameter D h = 4δ, wall temperatures and pressure. Cases without radiation A, B, C and D correspond to cases C1, C3, C4 and C2 in Ref. (Zhang et al. 2013) respectively, and radiative cases A_RAD, B_RAD, C_RAD, D_RAD03 and D_RAD08 to cases C1R1, C3R1, C4R1, C2R2 and C2R1 respectively.

  is the source term which acts as a pressure gradient term and drives the flow to obtain a target practical Reynolds number Re t D h . k is the turbulent kinetic energy and ǫ is the dissipation rate of k. The turbulent viscosity µ t is computed as C µ = 0.09, σ k = 1.4, σ ǫ = 1.4, C ǫ1 = 1.5, C ǫ2 = 1.9. The wall boundary conditions for Eq. (7.14), (7.15) and (7.16) are: u w = 0.0, k w = 0.0, ǫ w = 2 ν w (∂k 1/2 /∂y) 2 w . Using the turbulent viscosity µ t obtained from k-ǫ model, the temperature field is then calculated by solving the energy equation:

	✐													✐
	"thesis_yufang" -2013/11/27 -19:31 -page 180 -#180
	✐													✐
	180		7.2. NUMERICAL MODEL FOR FLOW FIELD
	y * =	ρ u ǫ y µ	with u ǫ = (νǫ) 1/4 ,			(7.19)
	Re t =	ρk 2 µǫ	.										(7.20)
	And the model function f ǫ writes
	f ǫ = 1 -exp -	y * 3.1		2	1 -0.3 exp -	Re t 6.5	2	(7.21)
	while the model constants are as follow: d dy c p µ Pr + µ t Pr t dT dy + S T = 0		(7.22)
	where the turbulent Prandtl number Pr t is modeled as (Kays 1994)
	is used, which resolves the following equations: Pr t = 1 0.5882 + 0.228(µ t /µ) -0.0441(µ t /µ) 2 (1 -exp( -5.165 µt/µ ))	.	(7.23)
	d dy and the prescribed wall temperatures are used as boundary conditions. A source (µ + µ t ) du dy + S u = 0, (7.14) term S T is added in Eq. (7.22) in order to sustain a target central temperature t value T c
	d dy	µ +	µ t σ k	dk dy	+ µ t	du dy	2	-ρǫ = 0,	(7.15)
	d dy	µ +	µ t σ ǫ	dǫ dy	+ C ǫ1	ǫ k	µ t	du dy	2	-C ǫ2 f ǫ	ρǫ 2 k	= 0,	(7.16)
	where S u µ t = C µ f µ	ρk 2 ǫ									(7.17)
	where
	f µ = 1 -exp -	y * 14	2	1 +	5 Re 3/4 t	exp -	Re t 200	2	,	(7.18)
	✐													✐
	✐													✐

  Taine. Direct numerical simulation of a turbulent channel flow coupled to radiative transfer. 23rd International Congress of Theoretical and Applied Mechanics, Beijing, China, 2012; • Y. F. Zhang, R. Vicquelin, O. Gicquel, J. Taine and Y. Huang. Wallmodeled LES with a new wall model accounting for radiation effects. 7th International Symposium on Radiative Transfer (RAD-13), Kusadasi, Turkey, 2013. "thesis_yufang" -2013/11/27 -19:31 -page 196 -#196 Appendix A -Analysis of the momentum equationThe equation of momentum in direction x, when averaged, becomes Moreover, for a periodic channel flow, one have: dρ u 2 dx = 0, ∂ρ u ′′2 ∂x = 0 and ∂ρ u ′′ w ′′ ∂z = 0 Hence the terms on the left-hand side remains only∂(ρ u ′′ v ′′ )Considering that the layer is thin, and the averaged values v and w are much smaller than u, the predominant term associated with the mean velocity is ∂ ∂y µ ∂u ∂y , while for the term with fluctuation, ∂ ∂y µ ′ ∂u ′ ∂y is the most important one. All the other terms are much smaller compared to these two terms. Moreover, as discussed in Ref.[START_REF] Cabrit | Direct simulations for wall modeling of multicomponent reacting compressible turbulent flows[END_REF]), ∂ ∂y µ ′ ∂u ′ Eq. (A.9) indicates that the flow inside a turbulent channel flow is driven by the gradient of pressure. In a periodic channel flow, since the pressure is constant, a source term is then required to play the same role as pressure gradient, i.e., to balance with the total stress gradient. τ tot has a linear profile when a spatially uniform source term is applied. Furthermore, in the close vicinity of walls, a constant total stress region can be assumed, indicating ′′ v ′′ ≈ 0.(A.11) This is the so-called equilibrium momentum equation for turbulent boundary layers.

												✐ ✐
	"thesis_yufang" -2013/11/27 -19:31 -page 197 -#197
												✐ ✐
	196	Appendix A -Derivation of turbulent boundary layer equations	197
	where the total stress τ tot writes
	∂(ρ u 2 ) ∂x τ tot = µ + ∂(ρ uv) ∂y du dy -ρ u ′′ v ′′ . + ∂(ρ uw) ∂z	= -	∂p ∂x	+	∂(τ 1j ) ∂x j	.	(A.5) (A.10)
	Since v = 0, w = 0, this equation can be reduced as
	∂[ρ( u 2 + u ′′2 )] ∂x	+	∂(ρ u ′′ v ′′ ) ∂y	+	∂(ρ u ′′ w ′′ ) ∂z	= -	∂p ∂x	+	∂τ 1j ∂x j	.	(A.6)
	∂y -ρ u The averaged form of the energy transport equation for a statistically stationary . The viscous shear stress term ∂τ 1j ∂x j writes ∂τ 1j ∂x j = ∂τ 11 ∂x + ∂τ 12 ∂y + ∂τ 13 ∂z = ∂ ∂x 2µ ∂u ∂x -2 3 µ ∂u ∂x + ∂v ∂y + ∂w ∂z (A.7) dτ tot dy = d dy µ du dy
	+ turbulent flow writes ∂ ∂y µ ∂u ∂y ∂(ρ uh ) ∂x + ∂(ρ vh) + ∂y + ∂v ∂x ∂(ρ wh) + ∂z ∂ ∂z = By separating the terms into the mean and fluctuation part, it becomes µ ∂u ∂z + ∂w . ∂x ∂ ∂x λ ∂T ∂x + ∂ ∂y λ ∂T ∂y + ∂ ∂z ∂τ 1j ∂x j = ∂ ∂x 2µ ∂u ∂x -2 3 µ ∂u ∂x + ∂v ∂y ∂w + u i ∂p ∂x i + τ ij ∂u i , ∂x j + ∂z	∂T ∂z (A.12) λ
	+ where the terms on the left hand can be separated into ∂ ∂y µ ∂u ∂y + µ ∂v ∂x + ∂ ∂z µ ∂u ∂z + µ ∂w ∂x + ∂ ∂x 2µ ′ ∂u ′ ∂x -2 3 µ ′ ∂u ′ ∂x + ∂v ′ ∂y + ∂w ′ ∂z ∂(ρ u h) ∂x + ∂(ρ u ′′ h ′′ ) ∂x + ∂(ρ v h) ∂y + ∂(ρ v ′′ h ′′ ) ∂y + ∂(ρ w h) ∂z	+	∂(ρ w ′′	(A.8)
			+	∂ ∂y		µ ′ ∂u ′ ∂y	+ µ ′ ∂v ′ ∂x		+	∂ ∂z	µ ′ ∂u ′ ∂z	+ µ ′ ∂w ′ ∂x	.
												∂y is negligible
	when compared to ∂ ∂y µ ∂u ∂y . Therefore, Eq. (A.6) is reduced to
	dp dx	=	d dy	µ	du dy	-ρ u ′′ v ′′ =	dτ tot dy	.	(A.9)
	✐ ✐											✐ ✐
	✐ ✐											✐ ✐

A.2 Analysis of the energy equation

  Appendix A -Analysis of the energy equationThis gives rise to the following simplified energy equation:For flows with a small Ma number (Ma< 0.2), the effects of pressure term and viscous term are negligible in the energy equation, hence, it yields

	✐												✐
	"thesis_yufang" -2013/11/27 -19:31 -page 198 -#198
	✐												✐
	198
	∂ ∂y	λ	∂T ∂y	-	∂(ρ v ′′ h ′′ ) ∂y	+ u i	∂p ∂x i	+ τ ij	∂u i ∂x j	= 0.	(A.15)
	dq tot dy	= -	d dy		λ	∂T ∂y	-ρ h ′′ s v ′′ = 0	(A.16)
	∂ ∂x		λ	∂T ∂x		+	∂ ∂y		λ	∂T ∂y	+	∂ ∂z	λ	∂T ∂z
	+	∂ ∂x	λ ′ ∂T ′ ∂x		+	∂ ∂y	λ ′ ∂T ′ ∂y	+	∂ ∂z	λ ′ ∂T ′ ∂z	(A.14)
	≈	∂ ∂y	λ	∂T ∂y		+	∂ ∂y	λ ′ ∂T ′ ∂y	≈	∂ ∂y	λ	∂T ∂y	.
	✐												✐
	✐												✐

  and p 2 are the uniform thermodynamic pressure and the hydrodynamic pressure respectively while p 2 /p 0 ≈ o(Ma 2 ).

	* 2 . The dimensional conservation equations for low-Mach number flows become ∂ρ ∂t + ∂(ρu i ) ∂x i = 0, (B.11) ∂(ρu i ) ∂t + ∂(ρu i u j ) ∂x j = -∂p 2 ∂x i + ∂τ ij ∂x j + S i i ∈ [1, 3], (B.12) ∂(ρh) ∂t + ∂(ρu j h) ∂x j = ∂p 0 ∂t -∂q j ∂x j + Q, (B.13) where p 0 "thesis_yufang" -2013/11/27 -19:31 -page 202 -#202

  τ ref , (D.2) t, n + 1 and n denote target values, prescribed time step n + 1 and n, respectively while subscript b is related to the bulk quantities as defined in section 4.3.1. The relaxation time τ ref is set to 0.3δ/u τ and V is the volume of the computational domain. The source term value is modified at each iteration when the simulation bulk Reynolds number is different from Re t b but tends to a constant value as a permanent regime is reached. However, in complex flows, the constant value that is obtained is different from S ref and the second term in Eq. D.2 remains constant and not null, necessarily introducing a finite bias between Re b and Re t b

	✐	✐
	"thesis_yufang" -2013/11/27 -19:31 -page 204 -#204
	✐	✐
	204
	where the superscript
	✐	✐
	✐	✐

  -ρ i ǫ i = S C k ,i + S P k ,i k i , (E.12)

				du	2
	,i			dy	i
	S ǫ,i = C ǫ1	ǫ i k i	µ t,i	du dy	2 i	-C ǫ2 f ǫ,i	ρ i ǫ 2 i k i	= S Cǫ,i + S Pǫ,i ǫ i ,	(E.13)
	with
	S C k ,i = µ t,i			du dy	2 i	; S P k ,i = -ρ i	ǫ i k i	;	(E.14)
	S Cǫ,i = C ǫ1	ǫ i k i	µ t,i		dy du		2

i ; S Pǫ,i = -C ǫ2 f ǫ,i ρ i ǫ i k i . (

E

.15)

If I - ν (yw, µ) were isotropic the first term of the second member of Eq. (6.22) would be similar to the classical result of the second one.

Remerciements

Acknowledgements

We thank the China Scholarship Council (CSC) for granting a doctoral fellowship to Yufang Zhang. This work was granted access to the HPC resources of [CCRT/CINES/IDRIS] under the allocations 2012-2b0164 and 2013-2b0164 made by GENCI (Grand Equipement National de Calcul Intensif).

Acknowledgements

We thank the China Scholarship Council (CSC) for granting a doctoral fellowship to the first author. This work was granted access to the HPC resources of [CCRT/CINES/IDRIS] under the allocations 2012-2b0164 and 2013-2b0164 made by GENCI (Grand Equipement National de Calcul Intensif). We are grateful to Dr. Philippe Riviere for providing the database of the Correlated-K model. We also acknowledge Dr. Vincent Moureau and the SUCCESS scientific group for providing the code YALES2 and support.

Abstract

In several conditions, radiation can modify the temperature law in turbulent boundary layers. In order to predict such an effect and the corresponding change in conductive heat flux at the wall, a new wall model for large-eddy simulation (LES) is proposed. The wall model describes the inner boundary layer which cannot be resolved by the LES. The radiative power source term is calculated from an analytical expression of the intensity field within the inner layer. In the outer layer, wall stress and conductive heat flux predicted by the wall model are fed back to the large-eddy simulation which is coupled to a reciprocal Monte-Carlo method to account for radiation. Several mixing-length models and turbulent Prandtl number formula are investigated. Then, the level of accuracy of the discretized radiation analytical model is investigated. Finally, fully coupled results are compared with Direct Numerical Simulation/Monte-Carlo results of turbulent channel flows at different Reynolds number, wall temperature and pressure conditions. The proposed wall model greatly improves the accuracy of the predicted temperature profiles and wall conductive heat fluxes compared to approaches without radiation accounted for in the inner layer.

Derivation of turbulent boundary layer equations

In this part, the thin equilibrium boundary layer equations are demonstrated for turbulent periodic channel flows. Both the momentum and energy equations are here considered. And these equations can also be derived similarly in most of the configurations.

A.1 Analysis of the momentum equation

For a turbulent flow in a statistically stationary state, the averaged equation of continuity (Eq. (3.2)) writes

Introducing the Favre average, it becomes

Since the periodic channel flow is homogeneous in z direction, it yields w = and ∂(ρ w) ∂z =. Moreover, when the flow is full developed, ∂(ρ u) ∂x = and hence, Eq. (A.2) becomes

An integration of this equation from the wall gives rise to

indicating that the Favre-averaged v is always zero in the whole channel flow.

If the density is variable, the Reynolds averaged wall-normal velocity v might not be zero [START_REF] Cabrit | Direct simulations for wall modeling of multicomponent reacting compressible turbulent flows[END_REF].

"thesis_yufang" -2013/11/27 -19:31 -page 201 -#201

Test of the adequacy of computational domain of channel flow

As discussed in section 3.2.2, in DNS, the computational domain should be large enough to represent the large turbulent scale properly. One way to ensure the adequacy of domain size in homogenize direction is to make sure that the fluctuations of main quantities are uncorrelated within half of the domain, i.e., the two-point correlation should decrease to zero within half size.

In present study, to illustrate the adequacy of the domain size (2πh, 2h, πh), the two point correlations for case C2 (described in section 4.3.5) are shown in Fig C .1 (the results in other cases are similar). It demonstrates that the domain is wide enough since most of the correlations decay to zero, except that the stramwise velocity correlation in the steamwise direction has a small value in the central region. In order to further examine the suitability of domain dimension in streamwise direction, a simulation of case C1 with a longer computational domain (4πh, 2h, πh) is performed and the results of u + and T + are compared in Fig C .2. The results of the two grids are the same (the results on the hot side are similar), hence, it is believed that the shorter domain size (2πh, 2h, πh) is appropriate and can be applied for all the cases. Let's assume the following definition for the source term S 1 :

where a 1 and a 2 are to be determined. Equation D.12 is then recast as

In this study, in order to obtain a relaxation towards the target value Re t b , a 2 is set as