Modelling turbulent heat flux accounting for Turbulence-Radiation Interactions

Abstract The present work investigates the modeling of turbulent heat transfer in flows where radiative and convective heat transfer are coupled. In high temperature radiatively participating flows, radiation is the most relevant heat transfer mechanism and, due to its non-locality, it causes counter intuitive interactions with the turbulent temperature field. These so-called Turbulence-Radiation Interactions (TRI) largely affect the temperature field, modifying substantially the turbulent heat transfer. Therefore, in the context of modeling (RANS/LES), these interactions require a closure model. This work provides the inclusion of TRI in the modeling of the turbulent heat transfer by adopting a unique approach which consists in approximating the fluctuations of the radiative field with temperature fluctuations only. Based on this approximation, coefficients of proportionality are employed in order to close the unknown terms in the relevant model equations. A closed form of all radiation-temperature-velocity correlation is explicitly derived depending on the chosen turbulent heat transfer model. This model is applied to a standard two-equation turbulent heat transfer closure and used to reproduce results obtained with high-fidelity DNS simulations. While a standard approach (i.e., neglecting TRI) is not able to correctly predict the DNS data, the new model’s results shows exceptional agreement with the high-fidelity data. This clearly proves the validity (and the necessity) of the proposed model in non-reactive, radiative turbulent flows.


Introduction
Many engineering applications work with high temperature fluids which are able to absorb and emit thermal radiation, such as H 2 O, CO 2 or CH 4 . As a consequence, the correct prediction of thermal radiative transport is of primary importance in high temperature application such as combustors, volumetric solar absorbers and heat transfer equipment in power plants. Radiative heat transfer is different from conduction due to its peculiarity of being inherently non-local. 2 This non locality causes counter intuitive interactions with other quantities such the as turbulent temperature field and conductive heat transfer.
All the studies regarding combined radiative and conductive heat transfer have reported the occurrence of turbulence radiation interactions (TRI) due to the highly non-linear coupling between temperature and the radiative transport. Namely, the turbulence behaviour of the temperature field causes the appearance of fluctuations of radiative quantities, which in turn modify the mixing in the flow and the relevant scales of the thermal structures.
The more practical works performed in the combustion field mainly deal with the effect of temperature fluctuations on the average radiative power and radiative transfer equation. The common conclusion is that this effect is negligible in non reactive flows, where the average radiative quantities can be calculated directly from the average temperature profile, but plays a large role in reactive flows where emission losses can be enhanced up to 30% due to TRI (Gupta et al., 2009(Gupta et al., , 2012Coelho, 2007Coelho, , 2012Coelho et al., 2003).
On the other hand, the first study which identified the role of radiation in the destruction of temperature fluctuations was performed by Townsend (1958). The study of the role of radiation in the modification of the temperature variance was advanced by findings in the field of atmospheric science (Schertzer and Simonin, 1982;Coantic and Simonin, 1984) which commonly deal with low temperature applications. In particular, Coantic and Simonin (1984) showed that radiative dissipation rate is proportional to κ p /ω K , where κ p is the Planck-mean absorption coefficient and ω K is the Kolmogorov wavenumber. Moreover, Soufiani (1991) has investigated the effect of radiation on the turbulent temperature spectrum of high temperature radiative gasses, demonstrating that radiative dissipation is more effective at the right end of the turbulent kinetic energy spectrum.
More recently, thanks to the increase in the available computational resources, several works have addressed TRI with direct numerical simulation (DNS) to investigate the structural change of temperature turbulence in presence of radiative heat transfer. Numerous works involve different approximations to solve the issue of high computational requirements for radiative calculation. The most used approximation is the optically thin approximation (OTA) which assumes an optically thin flow and, therefore, a negligible depletion of incoming intensity due to absorption. Employing this approximation Sakurai et al. (2010) investigated a horizontal buoyant radiative channel flow, discovering that radiative heat transfer leads to a destruction of the large organized buoyant structures. Other studies employed the gray gas approximation (Silvestri et al., 2018;Ghosh et al., 2014;Gupta et al., 2009) to visualize the effect of radiative heat transfer on temperature. In particular, our previous work (Silvestri et al., 2018) showed that TRI has a non linear dependency on optical thickness τ, which is caused by the contrasting roles of radiative emission and absorption.
Building on these results, we were recently able to parameterize TRI for a wide range of optical thicknesses in gray and non-gray participating turbulent flows . This allowed us to compare non-gray gasses in non-reactive flows to the much simpler gray gasses and accurately assess TRI directly using the values of the temperature variance and the information of the temperature length scales in the domain. Contrarily to previous works, we highlighted the dependency of radiative dissipation on the ratio κ g /ω c , where κ g is a TRIequivalent absorption coefficient and ω c is a characteristic wavenumber which accounts for anisotropic turbulent structures.
While all this knowledge has been gathered regarding the TRI mechanism and the coupling of convective and radiative heat transfer, standard models have not yet been adapted to these recent findings. When solving the Reynolds averaged Navier-Stokes (RANS) equations it is necessary to model the unclosed terms. It has been extensively proven that TRI modifies severely the value of the turbulent heat transfer, which is one of the variable requiring a proper closure model. This closure has never been provided in the presence of radiative heat transfer leading to the failure of all the developed models in high temperature, participating turbulent flows. Here we provide the inclusion of TRI in the modeling of the turbulent heat transfer following a rigorous mathematical procedure. We construct the model on the base of the knowledge gathered in our previous DNS investigations . This "TRI" closure model is applied to a standard two equation model and tested against various DNS cases to demonstrate both its necessity and validity.

Governing equations
As discussed later, the model is developed in the framework of a 1D turbulent channel flow enclosed between a hot and a cold wall due to the availability of high-fidelity DNS data. Nonetheless, the formulation is general and applicable to other flows and boundary conditions. The nondimensional Favre averaged Navier-Stokes equations for an emitting absorbing turbulent medium, read ∂ρ ∂t ∂ρθ ∂t where We have here assumed ideal gas, low Mach number, and constant heat capacity c * p . The asterisk * indicates a dimensional quantity, while the variables without asterisk are non-dimensional. In the above equations and in the rest of the text, the overbar and the tilde denote Reynolds and Favre averages, while the prime and double prime indicate their respective fluctuations. The non-dimensionalization is performed as following where x, t, u, ρ, μ, λ, θ, I and κ are non-dimensional position vector, time coordinate, velocity vector, density, viscosity, thermal conductivity, temperature, radiative intensity and absorption coefficient, respectively. Furthermore, U * b is the bulk velocity and quantities with subscripts c and h refer to values at the cold and hot wall, respectively. The nondimensional parameters appearing in Eqs. (1)-(3) are the Peclet number Pe = RePr and the Radiation number Rd = RePrPl, where Re is the bulk Reynolds number, Pr is the Prandtl number and Pl is the Planck number, defined as The radiative quantities of interest are absorption coefficient κ p , emission E and reference incident radiation G. If the spectral absorption coefficient κ * ν is not a function of the radiative wavenumber ν, the flow is considered to be gray and κ p = δ * P ( T * − 1 ) , where P represents a fifth order polynomial fit for the Planck-mean absorption coefficient of water vapour. Otherwise, κ p is calculated as in (Modest, 2013) Non dimensional emission and incident radiation are calculated as where T 0 is defined as . If the flow is gray, the standard definition of incident radiation is retrieved G = π − 1 ∫ 4π IdΩ. The Reynolds stress ρu ′′ i u ′′ j and the turbulent heat flux ρu ′′ j θ ′′ in Eqs. (2) and (3), respectively, are unclosed terms that require a closure. The modeling of the Reynolds stress is well established and, since radiative heat transfer does not directly affect velocity, it is not a topic of this study. On the other hand, the turbulent heat flux ρu ′′ j θ ′′ is greatly affected by radiative heat transfer and, therefore, its effect has to be accounted for to ensure correct temperature predictions.

Turbulence radiation interactions
Before deriving the complete model for the unclosed terms, we describe the most prominent interactions between radiation and turbulence in a participating media. Fig. 1 shows a schematic that describes all the interaction pathways between mean and fluctuating temperature and radiative quantities. Turbulence radiation interactions are divided into two main effects: the action of the radiation field (both mean and fluctuating) on temperature fluctuations (κ p , E, G→θ ′ ) and the development of a fluctuating radiative field by the action of temperature fluctuations (θ ′ →κ ′ p , E ′ , G ′ ). The two gray lines show pathways which are considered negligible in this study. The first one, which represents the modification of the mean radiative field due to the presence of temperature fluctuations (θ ′ →κ p , E, G), is relatively small in non reactive flows, where T ′ /T ≪ 1 (weak temperature fluctuations), as reported in several previous investigations (Gupta et al., 2009;Vicquelin et al., 2014;Roger et al., 2011;Ghosh et al., 2014;Silvestri et al., 2018). The second (κ ′ ,E ′ ,G ′ →θ) which quantifies the direct influence of radiative quantities fluctuations on average temperature, is much smaller than the counterpart (κ p , E, G→θ, average radiative field impact on mean temperature) and null in constant absorption coefficient flows . This pathway can be, therefore, safely neglected. The dominant effect of radiative fluctuations on θ is felt through the modification of θ ′ and the consequent change in turbulent heat transfer. This last pathway (κ p ,E,G ↔ θ ′ ↔ θ), which includes all the blue lines of Fig. 1, is considered here the most influential in non reactive flows and is hence the only focus of this study.
The specific novelty of this work consists in the modeling of all the relevant TRI pathways (blue lines) and the inclusion of these effects into a general turbulent heat transfer closure model (red lines). To avoid errors in modeling the radiative heat source, we will use the averaged radiative quantities (E, G, κ p ) directly from the reference DNS simulations. In order to model TRI, it is necessary to approximate the fluctuations of radiative quantities corresponding to the bottom right block in Fig. 1. For this purpose, we use the linear expressions recently derived in our previous DNS investigation  which relate the fluctuations of the radiative quantities (κ ′ , E ′ and G ′ ) to the temperature fluctuations θ ′ . Namely, where f κ , f E and f G are coefficients of proportionality, functions of the averaged quantities only, defined as following where ΔT * is T * h − T * c and c * 0 − c * 5 follow from the fitting of κ p . The TRIequivalent absorption coefficient κ g , appearing in the above expression, is defined as If non-gray, κ g in expression (10) depends on the averaged temperature profile as well as ω c which is a "characteristic" wavenumber that includes turbulence effects in the TRI model. In this case, κ g is found by averaging TRI effects over the whole absorption spectrum by considering each wavelength equal to an independent gray gas. Both ω c and κ g are derived in , but also discussed in more details in Section 5.2. For the full derivation of the expressions of f κ , f E , f G and κ g the reader is referred to our previous work .

Turbulent heat transfer two equation closure model
Most of the turbulent heat transfer closure models used in a RANS framework are based on the gradient-diffusion hypothesis which states that where α t is the flow dependent "eddy diffusivity". This quantity can be approximated in several ways. The most common is to relate it to the eddy viscosity μ t as α t = μ t /Pr t , with Pr t as the "turbulent Prandtl number", usually taken as constant equal to 0.9 (Antonia and Kim, 1991). This is a quite crude approximation which works in the limit of high Reynolds number flows when Pr ≈ 1. We have already demonstrated in our study on TRI in gray gases (Silvestri et al., 2018) that Pr t is largely modified by radiative heat transfer and cannot be used if the flow is able to emit and absorb radiation. The two equation model, on the other hand, does not rely on the turbulent Prandtl number, and estimates the turbulent diffusivity α t by relating it to a mixed time scale τ m , which incorporates both the velocity field and the temperature field information as where τ u is a time scale characteristic of the velocity field, while τ s represents the time scale of the thermal field. Usually, the contribution of these two time scales to the mixed time scale is considered equal (i.e., the exponents m and n are taken as m = n = 0.5) (Deng et al., 2001;Nagano and Kim, 1988;Sanders and Gokalp, 1998;So, 2000). In this work we followed the same approach. The expression for the eddy diffusivity then becomes where f m k is a damping function that accounts for low Reynolds number effects. This model has been developed and tested, with different details, in Deng et al. (2001), Nagano and Kim (1988), Sanders and Gokalp (1998) and So (2000). In this work we follow the model version developed by Deng et al. (2001). The model functions and constants involved are summarized later in Section 6. In order to assess the thermal time scale τ s , two additional non-dimensional transport equations, for temperature variance θ ′ 2 and dissipation of temperature variance ε θ , respectively, are solved (from here the name of the closure model). In contrast to previous works, we have here derived these equations for a radiatively participating flow, to account for the effect of radiative heat The blue lines show the pathways which will be accounted for in the present TRI model and, finally, the gray lines show negligible pathways in non-reactive flows (Roger et al., 2011;Silvestri et al., 2018;Vicquelin et al., 2014). On the other hand, dashed lines (which include the negligible terms) show pathways which are automatically accounted for since κ p , E and G are an input from the DNS solutions.
transfer. Both additional equations are displayed and discussed in the sections below.

Temperature variance transport equation
The exact transport equation for the temperature variance reads: In order to obtain a closed form of the above equation, several assumptions and approximations have to be applied. A first approximation consists in neglecting the terms containing θ ′′ (D θ in the above equation). This is exact in an incompressible flow, and, given the Morkovin hypothesis (Morkovin, 1962), which states that compressibility can be accounted for by considering mean density variations alone, it is an accurate approximation in low Mach number flows as well. This approximation implies that The validity of this assumption for the investigated cases is demonstrated by showing profiles of θ,θ, θ ′ 2 and θ ′′2 in Appendix A. In addition, it is common practice to assume that thermal conductivity fluctuations are low compared to its mean value (i.e., λ ≈ λ). Finally, the transport term is commonly modelled using a gradient-diffusion hypothesis employing the eddy diffusivity α t scaled by a coefficient σ θ . All these approximations lead to where P θ is the production of temperature variance, estimated, consistently with the turbulent heat flux, as and ε θ is the scalar dissipation calculated by its own transport equation.
The remaining unclosed term in this equation is the new radiative term R θ .

Scalar dissipation transport equation
The additional transport equation for ε θ has been derived for constant property flow and applied also to variable density flow keeping in mind that density fluctuations are low compared to the average density. For the case without radiation, the following formulation is the same as used in (Deng et al., 2001) ρ Again, the turbulent transport term has been modelled by employing the gradient-diffusion approximation as for the temperature variance transport equation. In Eq. (18) the terms P εθ and ε εθ are the production and the dissipation term, respectively. The first is modeled with the use of the mixed time scale, while the latter is divided into two different parts, one that accounts for velocity time scale and one which accounts for temperature time scale. Constants C p , C d1 and C d2 and model functions f p , f d1 , f d2 are taken as in Deng et al. (2001) and shown in Section 6. A new unclosed radiative term R εθ appears. The rationale for the derivation of this new term follows the procedure which leads to the ε θ budget equation. If we start from an incompressible, constant viscosity formulation, then In this case, ∂ t ε θ is equivalent to Therefore, the procedure to obtain the radiative term in the scalar dissipation budget equation follows these steps: (1) subtract the mean radiative heat source to the instantaneous radiative source, (2) derive the result in x j , (3) multiply by 2Pe − 1 ∂ xj θ ′ , (4) Reynolds average the resulting term. If stemming from the compressible low Mach number Navier-Stokes equations (expressed in a non-conservative form), ∂ t θ requires a division by ρ on the RHS. The first step of the procedure would then yield On the other hand, to account for variable density, it is common practice to multiply the scalar dissipation equation with ρ (Sanders and Gokalp, 1998). Therefore, by assuming that, in this context, the radiative terms obtained following the incompressible procedure is assumed to be valid also in a variable density framework. It has to be reminded that this assumption is valid only under the Morkovin hypothesis of weak density fluctuations and, therefore, it cannot be directly applied to high Mach number flows. To account for variable thermal conductivity, sinceassuming weak thermal conductivity final radiative term, as anticipated in Eq. (18) reads

TRI modeling
By following the approach presented in Silvestri et al. (2019), it is possible to find a closure for the additional radiative terms in Eqs. (16) and (18) by using the model functions shown in Eqs. (7)-(9). The radiative term in Eq. (16) can be rewritten in terms of emission and incident radiation by substituting Q = E − G. Additionally, a Reynolds decomposition of κ p , E and G yields From the modelling of radiative fluctuations as given in Eq. (16) Since f κ ≪κ p (absorption coefficient fluctuations are mostly negligible, see  and θ ′ 3 ≪θ ′ 2 , it is possible to safely neglect the last term on the RHS. The final model for the radiative term in the temperature variance budget equation reads The above equation is closed, as it depends on quantities readily available in a RANS framework (provided θ ′ 2 is modeled).
The radiative term R εθ in Eq. (18) can be expanded, by performing a Reynolds decomposition of κ p and Q, as As already done for Eq.
(3), we neglect the higher order term containing κ ′ p Q ′ . This is substantiated by the low impact of absorption coefficient fluctuations, demonstrated in . By further splitting the derivatives, yields the model for the radiative term in the ε θ transport equation Also the above expression, provided θ ′ 2 and ε θ are modeled, depends on quantities readily available and can be directly implemented in the turbulent heat transfer model. Note that κ p , E and G are taken from the correspective DNS simulations. However, our previous DNS investigation showed that the terms involving E and G in Eqs. (26) and (29) are negligible irrespective of the optical thickness value.

Modified temperature time scale
The definition of τ s to be used in the expression of the eddy diffusivity (13) accounts only for conductive dissipation of temperature variance (ε θ ). However, TRI acts as an additional "radiative dissipation" ε r which reduces drastically the temperature time scale (Townsend, 1958).
Therefore, to include TRI in the definition of α t , we must define a "modified" temperature time scale τ ★ s which accounts for ε r where C r is a model constant. Given the validation cases, the value of 0.5 for C r has been found to best match the DNS data. Radiative dissipation ε r , in a strict sense, is the dissipative part of the radiative term R θ which can be retrieved by expressing the radiative heat source in terms of divergence of radiative heat flux, By employing this definition it is possible to decompose the radiative term R θ into a dissipation ε r and a transport term ϕ r as On the other hand, we can assume that, away from the walls, the dissipation term is much larger than the transport term (ε r ≫ϕ r ) as shown in Silvestri et al. (2018). Therefore, it is possible to assume that ε r ≈ 0.5R θ . From this definition, the eddy diffusivity can then be corrected as follows It has to be pointed out that τ s is present also in the scalar dissipation transport Eq. (18) to model the dissipation and the production term (ε εθ and P εθ , respectively). Here, the original definition of τ s is maintained, as the influence of radiation is directly modeled through the term R εθ .

Characteristic wavenumber
The model employed for incident radiation fluctuations (G requires the estimation of a "characteristic wavenumber" ω c , which represents the length scale of the average energy-containing temperature structure. In  it is defined, such that anisotropy due to wall turbulence is accounted for, as where S θ is the one-dimensional temperature power spectrum. Since the power spectrum is not available in a RANS simulation, the integral length scale of temperature is used, calculated as in (Coantic and Simonin, 1984) In non-gray gas cases, κ g depends on ω c . An iterative procedure is, therefore necessary: κ g is initialized with κ p and, as the solution is iterated, κ g is updated as a function of the integral length scale and the mean dimensional temperature κ g l θ ⋅atan ( by integrating over line-by-line spectra retrieved from a high resolution, accurate spectral database (Rothman et al., 2012). For a detailed derivation of Eq. (36), the reader is referred to .

Summary of the model equations
In this section, the radiative modification is tested on different cases for which DNS data is available. The validation of the model is done by comparison with high fidelity DNS of a fully developed absorbing emitting turbulent channel flows enclosed by a hot and a cold wall which were performed in previous works (Silvestri et al., 2018. The DNS data is statistically steady and homogeneous in the streamwise and spanwise direction, and, as a consequence the model equations simplify to a one-dimensional problem (only the gradients in the wallnormal direction remain). Below, we summarize the model equations as well as the values for all the constants involved. It is reminded that we neglect fluctuation of transport properties (λ ≈ λ, μ ≈ μ) and make use of the gradient-diffusion hypothesis ρu ′′ v ′′ = − μ t ∂u ∂y and ρv ′′ θ ′′ = − α t ∂θ ∂y .
the model RANS equations for a statistically fully developed turbulent channel flow reduce to, ∂ ∂y ∂ ∂y Given the moderate Reynolds number of the test cases, turbulent vis-cosity is calculated using the v 2 − f model of Durbin (1995) which is able to correctly predict wall damping by introducing ad hoc damping relations. The model is not shown here, for more details the reader is referred to Durbin (1995). In particular, for variable density cases, the variable property formulation of Otero et al. (2018) is implemented, as it slightly improves the turbulent stress prediction (proof in Section 8.1).
The turbulent heat flux model equations are summarized below with the model functions and constants (as in Deng et al., 2001) shown in Table 1, and below

Test cases
The investigated cases are presented in Table 2. All cases are forced convection in a periodic channel bounded by an isothermal hot and cold wall. Both walls are black (∊ w = 1). In terms of thermal boundary conditions, the first 14 cases have T * h = 955 K and T * c = 573 K corresponding to T 0 = 1.5, while the last two have T * h = 1800 K and T * c = 600 K (T 0 = 0.5). All cases have a constant Planck number equal to 0.03. The DNS database includes constant and variable properties as well as gray and non-gray cases. In particular, the first seven cases are constant property, constant absorption coefficient and gray. They differ only in the magnitude of κ and, therefore, optical thickness (τ). These DNS cases are presented and discussed in our previous work (Silvestri et al., 2018)   where we characterized the effect of optical thickness on TRI.
The next four cases, designated by a ρ in the name, are still gray, but with temperature-dependent density and absorption coefficient. In particular the absorption coefficient is a 5 th order polinomial of T * − 1 , defined as where the constants c * 0 − c * 5 are taken from the Sandia national laboratories' (Barlow et al., 2001) model of water vapour Planck-mean absorption coefficient.
Finally, the last five cases have a spectrally varying, temperature dependent absorption coefficient and variable density. Fig. 2 shows the spectra for these cases at 800 K and 1 atm. These are (1) the spectra of water vapour named H 2 O, (2) the spectra of 10% carbon dioxide and 90% nitrogen, labelled CO 2 and (3) a synthetic spectra which mimics a multiphase medium, called QG for quasi-gray. The latter is generated by adding to the spectra of 10% H 2 O a constant (which represents the gray absorption coefficient of particulate media). The last two cases in Table 2 have a higher Reynolds number and variable viscosity and thermal conductivity.
Among all these cases, three transparent benchmarks (bench, benchρ and bench-highRe) are used to test the employed RANS models. The reference DNS results for the constant property cases have been obtained using the code described in Silvestri et al. (2018). The variable properties DNS results, on the other hand, are obtained by coupling the CLAM Finite Volume Method described in Silvestri et al. (2018) with the low Mach number Navier-Stokes solver described in Patel et al. (2015). Finally, the non-gray gas DNS results are produced by coupling the low Mach number Navier-Stokes solver with a high accuracy Monte Carlo spectral radiative solver described in .
In the rest of the section the following RANS model combinations will be compared: • v 2 − f for μ t with α t = μ t /0.9 • v 2 − f for μ t and θ ′ 2 − ε θ for α t with no TRI model (R θ = R εθ = 0) • v 2 − f for μ t and θ ′ 2 − ε θ for α t with TRI model We remind that, since the focus of this report is the modeling of turbulent heat transfer in presence of radiation, to avoid errors in the calculation of the radiative sources, the profiles of the average radiative quantities (κ p , E and G) are taken directly from DNS calculations. This ensures that, even if negligible, the θ ′ →κ p , E, G pathway is still accounted for (Fig. 1).

Transparent cases
The models are first tested on the transparent benchmarks to ensure correct implementation. Figs. 3 and 4 show calculated velocity and shear stress, and mean temperature and turbulent heat flux, respectively. Here, the dashed lines show the results obtained with a simple constant turbulent Prandtl number, while the solid lines are the calculations using the two equation turbulent heat flux model. In particular, the red lines show the results obtained using a classical implementation, based on Fig. 2. Spectra of the non gray cases at 800 K and 1 atm. H 2 O is the spectrum of water vapour, CO 2 is the spectrum of 10% carbon dioxide in N 2 , while QG is a synthetic spectra that mimics a multiphase medium. wall scaling, while the blue lines are obtained using an improved v 2 − f model, which accounts for variable properties, developed on the basis of semi-local scaling (Pecnik and Patel, 2017) and implemented in Otero et al. (2018). In the bench case, the turbulent heat flux model does not affect the velocity field, as temperature is a passive scalar. Since the results from all turbulent heat flux models would collapse on the same line, only one RANS result for u and ρu ′′ v ′′ is shown. It is possible to notice that the v 2 − f model slightly underpredicts the turbulent stress. For the bench-ρ and the bench-highRe cases, the classical and improved v 2 − f formulations for μ t are compared. It is possible to notice that, again, the choice of the turbulent heat flux model does not affect the velocity field (constant turbulent Prandtl number yields same results as the-two equation model), but the v 2 − f formulation does, with slightly improved results using the "semi-local scaling" implementation described in Otero et al. (2018) when large property variations are present (case bench-highRe). In the latter case, the turbulent stress is largely underpredicted on the hot side. The reason for this underprediction might be the very low density and high viscosity on the hot side that cause local low Reynolds number effects which are known to reduce the accuracy of a k − ε based turbulence model.
In constrast, the turbulent heat flux and mean temperature profiles (Fig. 4) show differences between the models used. If the constant turbulent Prandtl model is used, the turbulent heat flux is mispredicted in the center of the channel, leading to improved mixing and a mean temperature profile which is lower on the hot side and higher on the cold side of the channel. On the other hand, the two equation model leads to an overprediction of the turbulent heat flux in the core of the channel. This is caused by a slight overprediction of the turbulent heat transfer (derivative of ρv ′′ θ ′′ ) in the thermal conductive layer. Despite this overprediction, the important quantity is the derivative of the turbulent heat flux, which is better predicted with the two-equation model compared to the constant turbulent Prandtl number approach. This is proven by the better agreement of the mean temperature profile. Note that, as shown by Eq. (39), in a transparent case the temperature profile is completely defined by α t , and, therefore, reflects the real performance of a turbulence heat flux model. Therefore, the two-equation model leads to a smaller "thermal boundary layer" than the actual DNS data but an overall good performance in terms of average temperature profile in the core of the channel. Again, it is possible to notice that the variable properties v 2 − f model yields slightly improved results in case of large property variations (case bench-highRe). For this reason, the improved variable property v 2 − f formulation is used for all the following simulations.

Constant property, gray cases
Fig. 5 shows the results obtained for constant properties, low to intermediate optical thickness cases. As already demonstrated in several previous studies (Silvestri et al., 2018;Deshmukh et al., 2008;Vicquelin et al., 2014;Gupta et al., 2009;Zhang et al., 2013) if the channel is optically thin, the influence of TRI is negligible. Indeed, it is possible to notice that, since TRI is negligible, the differences between the models, for case gray-01, closely resemble the transparent cases. In particular, by assuming a constant turbulent Prandtl number, turbulent heat transfer is mispredicted in the center of the channel leading to a higher mean temperature on the cold side of the channel. The two equation model overpredicts the turbulent heat flux as in benchmark case, while the addition of a TRI model improves slightly the predictions. This does not translate in a visible improvement in the mean temperature profile as TRI impact is still very low. On the other hand, at an intermediate optical thickness TRI starts to play an important role, strongly affecting the turbulent temperature field. This influence is reflected in the failure of the standard models in predicting both the turbulent heat transfer and the average temperature field. In particular, turbulent heat flux is always severely over-predicted leading to an increased temperature mixing when compared to the DNS results. This is caused by the fact that standard models do not take into account the additional dissipative effect of radiative heat transfer on temperature fluctuations and, therefore, predict much higher thermal turbulence levels. In reality, TRI decouples the connection between turbulent velocity and temperature field, and severely reduce temperature fluctuations thanks to a longrange heat transfer mechanism. This results in a much lower turbulent heat flux compared to a non-radiative case with the same mean temperature gradient and velocity fluctuations. Contrarily, including a TRI model allows to predict the reduction in turbulent heat flux leading to accurate results in terms of mean temperature profile. It is possible to notice that increasing optical thickness (case gray-5) leads to a more severe misprediction of temperature by the standard models due to a higher TRI influence. Fig. 6 shows the optical thicker cases (gray-10, gray-20 and gray-10p). Again, since the impact of TRI on the turbulent temperature field is very high, the use of the proposed model is necessary to achieve an accurate prediction. Nonetheless, we imagine that substantially increasing the optical thickness (τ≫1) would decrease the necessity of the TRI model. This is caused by the fact that the strength of TRI is not a monotonic function of τ, it is zero at τ = 0, reaches a maximum and returns to zero at τ→∞. Case gray-10p shows the results obtained for a lower Prandtl number (Pr = 0.7). The TRI model performs still  exceptionally. We did not test Prandtl numbers larger than unity because of the unavailability of DNS data. Nevertheless, since the relative influence of the radiative heat transfer (and TRI) is reduced with an increase in Prandtl number, we believe that the results will still be accurate without any modification. It is possible to imagine that for a higher Prandtl number (Pr > 1) temperature structures will be significantly different, and this change will have to be taken into account with a Prandtl number dependency of the characteristic wavenumber ω c (approximated by l θ ). Nonetheless, the strength of TRI, which scales with Pr − 1 , will reduce in intensity. Therefore, it is straight-forward to show that also for an increased Prandtl number, no significant modification of the proposed model would be required.
Finally, it is important to notice how, in both Figs. 5 and 6, the two equation model without TRI performs very poorly in terms of turbulent heat transfer. We already showed that the two equation model results in an over-estimation of turbulent heat transfer near the walls (Fig. 3). In case of radiative heat transfer, where turbulent heat transfer is largely suppressed near the boundaries, this over-estimation becomes unacceptable. TRI, especially in case of an intermediate to high optical thickness, is the dominant mechanism. Thus, by including a closure model for TRI we are able to correct the mispredictions and obtain excellent results. The constant turbulent Prandtl number model (red dashed lines) did not necessarily over-predict the turbulent heat flux near the walls in the benchmark cases, but is still a very crude approximation which connects tightly the velocity and the temperature field. This connection is partly severed by TRI leading also to unacceptable results as optical thickness (and radiative heat transfer strength) increases. Fig. 7 shows results for the variable property, gray cases compared to DNS. As explained in  the fluctuation of absorption coefficient do not impact TRI significantly as much as E and G fluctuations. As a consequence, the results for these cases follow very closely what was observed in Figs. 5 and 6. In particular, the low optical thickness case does not feel the impact of a fluctuating radiative field  and adding a TRI closure model does not change significantly the results. Contrarily, a large improvement is found when including TRI for intermediate and large optical thickness cases (gray-ρ1 and gray-ρ10). To show the accuracy of the length scale approximation, Fig. 8 shows profiles of ω c as calculated from DNS data with Eq. (34) and l − 1 θ approximated with Eq. (35) in the RANS simulations. It can be noticed that near the walls l − 1 θ has an unphysical spike which is not seen in the DNS data. While ω c has a fixed value at a boundary, l θ tends to zero as k goes to zero, which means that, following expression (35), the integral length scale of temperature structures tends to zero as a fixed boundary is approached. This caused the characteristic wavenumber calculated by the RANS model to tend to infinity. Another misprediction is seen in the value of ω c for the optically thickest case gray-ρ10. This is attributed to the fact that Eq. (35) is not able to predict thermal structure's enlargement caused by the dual absorption-emission process at high absorption rates (Silvestri et al., 2018) as it does not include any radiative quantity. Nonetheless, both of these deviations from DNS values (which will be addressed later in Section 8.5 while discussing second order statistics) do not impact the predictions of both turbulent heat transfer and mean temperature.

Non-gray cases
Figs. 9 and 10 present the results for the non-gray cases with a spectrally varying absorption coefficient. Fig. 9 shows, from left to right, cases with an H 2 O, CO 2 and QG type spectrum (see Fig. 2). On the other hand, results presented in Fig. 10 are obtained with an H 2 O type spectrum. For the non-gray cases, the failure of the standard models is not as trivially connected to the optical thickness shown in Table 2 as for the previously analysed gray cases (i.e., for gray cases higher τ lead to larger misprediction if the TRI model was not included). In particular, case spec-Part (with τ p = 2.79) shows the worst performance in terms of both turbulent heat flux and mean temperature. On the other hand, case spec-    CO 2 , with a comparable optical thickness (2.99) seems to be relatively unaffected by TRI. This is because, as we explain in , two different parameters control TRI in a non-gray gas: τ p which is the Planck-mean optical thickness, and τ g , which is the TRI-equivalent optical thickness, based on the parameter κ g . These two parameters are defined as follows τ p = 0.5 ∫ 2 0 κ p dy, and τ g = 0.5 In gray gas cases τ = τ g = τ p , meaning that TRI depends only on the Planck-mean optical thickness. This is not the case in non-gray gases. The two optical thicknesses are shown, for the different cases, in the Table 3, together with the gray-equivalent optical thickness (τ eq ) described below. Since the influence of absorption coefficient fluctuations is generally negligible , the impact of TRI on the temperature field scales with where κ g ⋅f ( κ g ) is bounded between 0 and 1 and increases with κ g . The RHS of the proportionality relation has the dimensions of an absorption coefficient and is obtained by dividing R θ by . Therefore, the impact of radiation increases with τ p and decreases with τ g . To compare the results shown in Fig. 9 with the gray gas cases it is possible to define a new optical thickness (we call it here gray-equivalent), which is a combination of these two parameters. If we assume that only one absorption coefficient (here denoted κ eq ) is responsible for TRI (as in gray gasses), the "gray-equivalent" optical thickness can be defined as: The values of τ eq for each case, calculated with DNS data and ω c as in Eq.
(34), are given in Table 3. Given these values, it is possible to compare the predictions obtained for the non-gray cases with the gray gas cases. For spec-H 2 O, (τ eq = 0.74) neglecting TRI results in slightly better prediction than for the gray cases with τ = 1 (gray-1 and gray-ρ1), but still unacceptable if compared to the optically thin gray cases (gray-01 and gray-ρ01). In this case, the inclusion of a TRI model is necessary to obtain satisfactory results. Moreover, spec-CO 2 , despite the seemingly high τ p , is very similar to the optically thin cases (gray-01 and gray-ρ01) given that τ eq is equal to 0.118. Therefore, the two equation model is improved only slightly when including the TRI closure. On the other hand, since spec-Part has very similar τ p and τ g (the absorption spectrum has a really low variability), the result is similar to an optically intermediate gray case with optical thickness between 1 and 5 and hence largely impacted by TRI. Finally, case highRe-H 2 O has a τ eq lower than case spec-H 2 O, but is more affected by TRI (i.e., the deviations of the standard models are larger). This is caused by the fact that, in this case, there is more emitted energy that can dissipate thermal fluctuations. In simple terms, the average dimensional temperature is higher or, mathematically, Nevertheless, the TRI closure model seems to yield remarkable results independently from the Planck-mean, TRI-equivalent or gray-equivalent optical thickness of the case. Fig. 11 shows the comparison between the integral thermal length scales (top) and the TRI-equivalent absorption coefficient (bottom), calculated with DNS data and obtained with the RANS simulations (the latter is shown in Fig. 10 for the highRe-H 2 O case). The red dashed lines are the Planck-mean absorption coefficient profiles, plotted for reference. As seen previously for the variable property cases, the integral length scale l θ tends to zero while approaching the boundary, causing a spike in the characteristic wavenumber l − 1 θ . This is reflected in the calculation of the TRI-equivalent absorption coefficient which also spikes near the walls. Fortunately, the relevant function which approximates the incident radiation fluctuation is dependent on the ratio κ g /ω c . Therefore, the misprediction of ω c is corrected by the calculation of the TRI-equivalent absorption coefficient. Aside the problems close to the boundaries, the model (and iterative approach) employed seems to yield fairly correct values, especially for the calculation of κ g , which, due to its dependency on both temperature and thermal length scales, is the most complex quantity to assess.

Second order statistics
In this section we present the quantities calculated by the two equation model (θ ′ 2 and ε θ ) and the radiative dissipation assessed by the additional TRI closure model (R θ ). Since, in case of an eddy diffusivity calculated using a constant turbulent Prandtl number these quantities are not available, only the two equation model is shown in comparison to DNS data. Fig. 12 shows the results obtained for the constant property gray cases. It is possible to notice that for a low to intermediate optical thickness (gray-01 and gray-1), the TRI model predicts very accurately radiative dissipation, leading to a satisfactory calculation of temperature variance and molecular dissipation. Already for τ = 1, not accounting for TRI causes temperature variance and molecular dissipation to be unphysically high. This is caused by the large mean temperature gradient which develops in the center of the channel that translates in a high temperature variance turbulent production rate. The high production rate, coupled to absence of a radiative dissipation model results in a large overprediction of temperature variance. Increasing the optical thickness further (gray-10) leads to a lower accuracy in terms of radiative dissipation. This is most likely caused by the overprediction of the characteristic wavenumber. Since, for the high optical thickness cases, the characteristic wavenumber is effectively higher than the actual value obtained from the DNS simulations (Fig. 8), predicted incident radiation fluctuations (G ′ ) are lower, see Eq. (9), which result in a higher absolute value of radiative dissipation, see Eq. (25). Physically speaking, the model is predicting smaller thermal structures which are optically thinner and not capable of re-absorbing emitted radiation. This will lead to a decreased temperature variance (and ultimately a lower equilibrium R θ ). Therefore, the cause for having a lower R θ in Fig. 12 (for high optical thickness cases) is, counter-intuitively, an overprediction of radiative dissipation. Nevertheless, the temperature variance obtained without accounting for TRI is progressively worse as the optical thickness increases due to the larger mean temperature gradient. Therefore, despite slight inaccuracies in the second order statistics, including TRI leads to excellent results also for higher optical thickness cases. Figs. 13 and 14 show the performance of the RANS simulations in terms of second order statistics for variable properties and non-gray cases, respectively. As for the average profiles, the conclusions drawn for constant property cases are suitable for the variable property gray cases with the same optical thickness. It is interesting to notice the high accuracy of the TRI model in an intermediate optical thickness scenario, as for gray-ρ1 and all the non-gray cases. Also for highRe-H 2 O, which has variable transport properties and a significantly larger ΔT * and Re, the TRI model approximates very accurately radiative dissipation, allowing a correct prediction of temperature variance and molecular dissipation.

Conclusions
We developed a general radiative modification which can be applied to most turbulent heat transfer models. The modification consists of a first order approximation of the fluctuating radiative field which is expressed as a linear function of temperature fluctuations. This TRI closure model is then applied to a two-equation turbulent heat flux model which evaluates temperature variance and scalar dissipation rate.
The improved model has been tested on several cases in comparison with available DNS data to prove its validity. The results show that in case of a radiative flow, the proposed model is always capable to improve the results when compared to the standard models available.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Verification of the Morkovin hypothesis
To assess the impact of density fluctuation on average quantities and second order statistics, plots of Reynolds and Favre averaged temperature and temperature variance, obtained by direct numerical simulation, are shown in Fig. A.15. It is clearly demonstrated that in these low Mach number cases, the Morkovin hypothesis is fulfilled.