Cédric Flageul
email: cedric.flageul@ijs.si

Iztok Tiselj
email: iztok.tiselj@ijs.si

Sofiane Benhamadouche
email: sofiane.benhamadouche@edf.fr

Martin Ferrand
email: martin.ferrand@edf.fr

Wall-resolved LES of channel flow: A correlation for the discontinuity of ε θ at the fluid-solid interface

Wall-resolved LES of channel flow: A correlation for the discontinuity of ϵθ at the fluid-solid interface

Introduction

Conjugate heat transfer refers to the thermal coupling between a fluid and a surrounding solid. It is of prime importance in industrial applications where thermal fatigue is a concern. In the nuclear field, thermal fatigue and fluctuating thermal stresses are particularly important in case of a severe emergency core cooling or long-term ageing of materials.

Such complex applications are often studied experimentally and numerically. However, numerical investigations of turbulent flows at very high Reynolds numbers remains challenging. CFD (Computational Fluid Dynamics) analysis of industrial applications usually relies on high-Reynolds RANS (Reynolds-Averaged Navier-Stokes) and occasionally on wall-modelled LES (Large Eddy Simulation).

The turbulent flow applies on the solid domain a thermal load characterized by a broad spectrum. Qualitatively, the heat diffusion in the solid domain provides a strong damping at higher frequencies. Therefore, smaller scales tend to apply a thermal stress in the vicinity of the fluid-solid interface while larger scales tend to penetrate deeper in the solid. Refined analysis actually shows that high stress amplitude events are generally associated with low probability events (Costa Garrido et al. [START_REF] Garrido | Uncertainties in the thermal fatigue assessment of pipes under turbulent fluid mixing using an improved spectral loading approach[END_REF]), thus making accurate estimation of thermal fatigue in industrial applications even more challenging as the CFD simulations should provide at least a few minutes of operation in realistic conditions.

Analytical studies on conjugate heat transfer in turbulent flows were pioneered by Polyakov [START_REF] Poliakov | Wall effect on temperature fluctuations in the viscous sublayer[END_REF] and Geshev [START_REF] Geshev | Influence of heat conduction of the wall on the turbulent prandtl number in the viscous sublayer[END_REF]. The fundamental solutions of the heat equation in the solid domain have a non-compact support. For instance, semi-infinite solids with a flat fluid-solid interface subjected to a statistically steady forcing can be characterized by a compatibility condition at the fluidsolid interface expressed as a spatio-temporal convolution (Flageul et al. [4]). Such non-local effects are specific to conjugate heat transfer and tend to become negligible only when the thermal properties of the fluid and of the solid differ by orders of magnitude.

Numerical studies on conjugate heat transfer in turbulent flows were pioneered by Kasagi et al. [START_REF] Kasagi | Numerical investigation of near-wall turbulent heat transfer taking into account the unsteady heat conduction in the solid wall[END_REF] and their 2D synthetic turbulence model. The first DNS (Direct Numerical Simulation) with conjugate heat transfer was a turbulent channel flow, performed by Tiselj et al. [START_REF] Tiselj | DNS of turbulent heat transfer in channel flow with heat conduction in the solid wall[END_REF]. Following those studies, some of the present authors and co-workers performed additional DNS of the turbulent channel flow [START_REF] Flageul | DNS of turbulent channel flow with conjugate heat transfer: Effect of thermal boundary conditions on the second moments and budgets[END_REF] to extract the budgets of the second moments (i.e. the turbulent heat fluxes and the temperature variance). This work was motivated by the global lack of validation data for second-order RANS turbulence models (Dehoux et al. [START_REF] Dehoux | An elliptic blending differential flux model for natural, mixed and forced convection[END_REF], [START_REF] Dehoux | Algebraic modeling of the turbulent heat fluxes using the elliptic blending approachapplication to forced and mixed convection regimes[END_REF]) in case of an imposed heat flux, in case of a heat exchange coefficient, or in case of conjugate heat transfer: most of the previous DNS of the turbulent channel flow have an imposed temperature at the wall.

To the best of the authors knowledge, the only RANS turbulence model designed to take into account conjugate heat transfer -i.e. able to solve the temperature fluctuations both in the fluid and in the solid -was published by Craft et al. [START_REF] Craft | Towards the development of RANS models for conjugate heat transfer[END_REF]. It is also the only RANS turbulence model designed to correctly handle cases with an imposed heat flux (Mangeon et al. [START_REF] Mangeon | Unifying the near-wall treatment of the turbulent heat fluxes for all kind of temperature boundary condition with the elliptic blending approach[END_REF]). However, it was recently shown that the dissipation rate (ε θ ) associated with the halved temperature variance (T 2 /2) is discontinuous at the fluid-solid interface in case of conjugate heat transfer (Flageul et al. [START_REF] Flageul | On the discontinuity of the dissipation rate associated with the temperature variance at the fluid-solid interface for cases with conjugate heat transfer[END_REF]). Although there is currently no coupled RANS model taking this discontinuity in account there is a global agreement that one is needed (Shams et al. [START_REF] Shams | Synthesis of a CFD benchmarking exercise for a T-junction with wall[END_REF]).

The background on such turbulence models is given in the paper of Craft et al. [START_REF] Craft | Towards the development of RANS models for conjugate heat transfer[END_REF] and in classic textbooks on turbulence modelling. A simple -and thus approximate -sketch would be to say that RANS models use the turbulent kinetic energy k and the associated dissipation rate ε to model velocity fluctuations. Regarding temperature fluctuations, the halved temperature variance T 2 /2 and the associated dissipation rate ε θ are very similar to k and ε, respectively. However, velocity fluctuations remain inside the fluid domain while temperature fluctuations may penetrate inside solid domains adjacent to the fluid one. In case of conjugate heat transfer, this fundamental difference must be taken into account, and the turbulence model should be able to evaluate T 2 and ε θ in the solid domain adjacent to the fluid one.

In a previous paper, we have demonstrated the ability of wall-resolved LES to estimate the discontinuity of ε θ on heated channel flows using Code Saturne, Électricité de France in-house and open-source CFD software (Flageul et al., [START_REF] Flageul | Quantification of the discontinuity of the temperature variance dissipation rate: wall-resolved LES of turbulent channel flow with conjugate heat transfer[END_REF]). In the present paper, we present deeper analysis in the form of a power-law correlation allowing one to predict this discontinuity of ε θ on channel flows and thus impose the appropriate treatment at the interface. Wall-resolved LES at various Prandtl and Reynolds numbers have been performed to assess the robustness of the proposed correlation. This is a step forward towards a rich validation database for future RANS models adapted to conjugate heat transfer (Mangeon et al.,[START_REF] Mangeon | Unifying the near-wall treatment of the turbulent heat fluxes for all kind of temperature boundary condition with the elliptic blending approach[END_REF]).

Governing Equations and Cases Description

Governing equations

We consider the turbulent flow of a Newtonian fluid with constant physical properties. Furthermore, we assume that the flow is incompressible and that the physical properties in the solid domain are constant. The subscripts f and s are used for the fluid and for the solid, respectively.

In this subsection, we omit boundary conditions for the sake of simplicity. Firstly, the conservation of mass in a fluid with a constant density reduces to the incompressibility condition

∂ i u i = 0 (1)
where u i is the velocity in direction i.

Secondly, the conservation of momentum in the fluid domain is

∂ t u i + ∂ j (u i u j ) = - 1 ρ f ∂ i P + ∂ j σ ij + f i ( 2
)
where ρ is the density, P the pressure, σ ij the deviatoric part of the stress tensor and f i a source term. We use the WALE subgrid-scale model (Nicoud and Ducros [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF]) as implemented in the version 5.0 of Code Saturne: (3

σ ij = (ν + ν t ) S ij with ν t = 4 √ 2V 2/3 cell C 2 w S d kl S d
)
where ν is the kinematic viscosity, V cell the volume of the computational cell, C w the model constant, S d ij the deviatoric part of the square of ∂ j u i and S ij the deviatoric part of ∂ j u i . The constant C w is set to 0.25, the default value in Code Saturne.

Lastly, the energy conservation equation is

∂ t T f + ∂ i (T f u i ) = ∂ j α f + ν t P r t ∂ j T f + f T f (4)
in the fluid domain,

∂ t T s = α s ∂ jj T s + f Ts (5)
in the solid domain and

T f = T s and λ f ∂ n T f = λ s ∂ n T s (6)
at the fluid-solid interface. Here, T is the temperature, α the thermal diffusivity, P r t the turbulent Prandtl number, f T a source term and λ the thermal conductivity. The turbulent Prandtl number is set to 0.5 in the present study (Grötzbach [15]).

As we investigate forced convection, the transported scalars -like the temperature -do not impact the velocity field or the pressure. Thus, we can transport simultaneously an arbitrary number of passive scalars, and each scalar can use different thermal properties in the fluid and in the solid. This allows significant savings in terms of CPU time: up to 51 simultaneous scalars are considered in our simulations. In addition, each passive scalar is transported by the same velocity field. This peculiarity has eased our investigation as it reduced the impact of the statistical uncertainty (Flageul and Tiselj [16]).

Geometry

Figure 1 is a sketch of the domain. x, y and z are the streamwise, wall-normal and spanwise directions, respectively. The fluid domain is bounded (-δ < y < δ) and δ is the channel halfheight. The solid domains on top and bottom of the fluid one are located at y > δ and y < -δ, respectively. Both have the same total height, δ S . The exact extension of the channel in the streamwise and spanwise directions and the height of the solid domains depends on the case considered, as described in the subsection 2.5. . Based on previous studies, the present authors estimate that for all the cases considered here and summarized in the Table 1, the solid domains extension in the wall-normal direction is sufficient so that the boundary condition used at the outer wall has no significant impact on the statistics at the fluid-solid interface.

For all the LES considered here, the size of the cells in wall-units in the streamwise and spanwise directions is δx + = 30 and δz + = 15, respectively. On both sides of the fluid-solid interface, the first cell has an extension in the wall-normal direction of 1 wall-unit. Further away from the wall, the wall-normal extension of the cell is bigger, according to a geometric law with a progression factor of 1.09. In the fluid, as soon as the wall-normal extension of the cell reaches 15 wall-units, this progression factor is set to 1 so that the wall-normal extension of the cells does not exceed their spanwise extension. All the LES use a constant time step with a maximal instantaneous CFL number below 0.5.

Source terms

The source term in the momentum equation is an imposed pressure gradient, which is constant in space and in time. It compensates the viscous friction at the wall and keeps the flow in a statistically steady state. Using the averaged momentum equation, the no-slip boundary condition and the vanishing ν t at the wall, one can derive δf x = u 2 τ where u τ is the friction velocity associated with the wall shear stress. The Reynolds number based on the friction velocity and the channel half-height is Re τ = uτ δ ν . This definition, combined with the previous relation, allows us to choose the source term f x so that the target Re τ is exactly reached.

The source term in the energy equation in the fluid domain depends on the instantaneous streamwise velocity and on the instantaneous bulk velocity, as defined by Kasagi et al. [START_REF] Kasagi | Direct numerical simulation of passive scalar field in a turbulent channel flow[END_REF]. It is a volumetric sink term which compensates the heat input at the boundaries, so that the case remains statistically steady. Using the previous notations and writing q w the heat flux per unit surface inputted at the top and bottom boundaries, the source terms are:

f T f = 2α f q w L x L z λ f Lx x=0 δ y=-δ Lz z=0 u x u x and f Ts = 0 (7)
The subsection 2.4.further describes the boundary conditions used.

The source terms allow one to derive a theoretical friction velocity and friction temperatures which are not plagued by any statistical uncertainty. In the present paper, those theoretical friction values are used to express raw statistics in wall-units.

At this stage, in addition to the friction Reynolds number, we can introduce additional dimensionless numbers. The Prandtl number is P r = ν α f . The fluid-to-solid thermal diffusivity ratio is G = α f αs . The solid-to-fluid thermal conductivity ratio is G 2 = λs λ f . The thermal activity ratio -which is also the fluid-to-solid thermal effusivity ratio -is

K = 1 G 2 √
G . Although we have 3 dimensionless numbers describing the fluid and solid thermal properties ratios, only 2 are independent (Tiselj and Cizelj [START_REF] Tiselj | DNS of turbulent channel flow with conjugate heat transfer at prandtl number 0.01[END_REF]).

Boundary conditions

For the streamwise and spanwise directions, periodicity is used. In the wall-normal direction, the 3 components of the velocity vanish at the fluid-solid interfaces. Regarding the transported scalars, different cases are envisaged. If the scalar is not coupled, either its value (Dirichlet boundary condition) or its flux (Neumann boundary condition) is imposed at the fluid boundaries, located at y = ±δ. Actually, an imposed value at the fluid boundary is usually a good approximation when simulating the temperature inside a fluid in contact with a conductive solid. An imposed flux would rather correspond to the temperature inside a fluid in contact with a solid of low thermal conductivity.

In case of an imposed value, the temperature is arbitrarily set to zero at the boundary. Such cases will be referred to as isoT . In case of an imposed flux, its value must match the volumetric heat sink term imposed in the fluid domain. Such cases will be referred to as isoQ. Obviously, if a scalar is not coupled, it does not correspond to conjugate heat transfer.

For a coupled scalar (i.e. in case of conjugate heat transfer), the flux is imposed at the outer walls, located at y = δ + δ S and y = -δδ S , respectively. Here again, the value of the flux must match the volumetric heat sink term imposed in the fluid domain. For coupled scalars, there is continuity of the scalar and of its diffusive flux at the fluid-solid interface. Regardless of the case, the same boundary condition is imposed on both sides, so symmetry is preserved.

Simulations summary

All the LES presented include 2 non-coupled scalars. Those non-coupled scalars actually represent the limit of coupled cases when the solid thermal conductivity is very high, or very low, compared to the fluid one. The Dirichlet and Neumann boundary conditions correspond to the cases λ s λ f and λ s λ f , respectively. All the simulations also include at least one coupled scalar, with the same thermal properties in the fluid and in the solid.

In the section 3. , four LES are presented. The Reynolds numbers investigated are 395 and 1020. The Prandtl numbers investigated are 0.71 and 1. Table 1 summarizes the performed simulations. The case Re τ = 395 and P r = 0.71 includes 49 coupled scalars: 7 values of the ratios G and K are simultaneously simulated. The other cases include 3 coupled scalars: one with unit fluid-solid thermal properties ratios, one with G = 1.3 and K = 2.8 and one with G = 0.1 and K = 0.23. The second case is an approximation of the thermal coupling of water and plexiglas. The last one is an approximation of the thermal coupling of water and steel. This is further discussed in section 3. .

Discretization

Code Saturne 1 is an open-source CFD solver for incompressible or weakly dilatable flows. It is based on the FVM library and can handle unstructured meshes. The finite volume solver is collocated. The predictor/corrector algorithm used for pressure-velocity coupling is combined with a Rhie and Chow filter to avoid odd-even oscillations. Further details can be found in Archambeau et al. [START_REF] Archambeau | Code saturne: A finite volume code for the computation of turbulent incompressible flows-industrial applications[END_REF].

For the velocity and the scalars, the convection scheme used in this study is fully centered: the existing slope-test by default for the scalars is deactivated. Regarding time advancement, a Crank-Nicolson scheme is used, except for the convective term which uses an Adams-Bashforth time scheme for the transporting velocity. As our meshes are made of orthogonal hexahedra, no gradient reconstruction sweeps are needed.

Extrapolation of statistics at the fluid-solid interface

First, we recall the definition of the dissipation rate:

ε θ = α∂ i T ∂ i T (8)
As shown in [START_REF] Flageul | On the discontinuity of the dissipation rate associated with the temperature variance at the fluid-solid interface for cases with conjugate heat transfer[END_REF], this quantity is discontinuous at the fluid-solid interface, where the ratio of the solid and fluid dissipation rates verify

ε θ,s ε θ,f = ∂ y T f ∂ y T f ∂ i T f ∂ i T f K 2 + 1 - ∂ y T f ∂ y T f ∂ i T f ∂ i T f 1 G (9)
In Code Saturne, most of the relevant physical quantities are defined at the center of the cells. However, the fluid-solid interface is located at the faces of the cells and not at their center. In this subsection, we describe the strategy used to extrapolate statistical quantities at the fluidsolid interface for coupled fields. This very important post-processing step is performed after the simulation and uses quantities averaged in time and over homogeneous directions. In the following, y f and y s are the coordinate of the center of the first fluid and solid cells, respectively, and y f s is the coordinate of the fluid-solid interface.

For the streamwise contribution to ε θ , one can use a first-order Taylor expansion. Combined with the continuity of temperature and heat flux at the fluid-solid interface, one gets

∂ x T f ∂ x T f (y) = ∂ x T f ∂ x T f (y f ) + b f (y -y f ) ∂ x T s ∂ x T s (y) = ∂ x T s ∂ x T s (y s ) + b s (y -y s ) ∂ x T f ∂ x T f (y f s ) = ∂ x T s ∂ x T s (y f s ) b f = G 2 b s ( 10
)
The situation is exactly the same for the spanwise contribution to ε θ . For the wall-normal contribution to ε θ , the situation is different as there is no compatibility condition for its derivative. Thus, we use the continuity at the fluid-solid interface of the onepoint correlation between the temperature and its wall-normal derivative. We define the angle φ with

cos (φ) = T ∂ y T T 2 ∂ y T ∂ y T (11)
The numerator and the first term in the denominator can be estimated with T 2 . In this study, cos (φ) at the fluid-solide interface is approximated with its value at y s . Although one can estimate it with a combination of the values at y f and y s , the authors have observed that such combinations tend to produce non-physical ε θ at the fluid-solid interface. For instance, estimating cos (φ) at y f s with the halved sum of the values at y f and y s can produce a value of ε θ,s lower at y f s than at y s . Such a result is not admissible: for a statistically steady case, without source term in the solid domain, the diffusion of ε θ,s must remain positive in the solid domain. For the present case, it implies a decaying ε θ,s with increasing distance to the fluid-solid interface.

To obtain the variance of the temperature (and its derivative) in the cells adjacent to the fluidsolid interface, one can use a second-order Taylor expansion. As the LES is wall-resolved and the flow statistically steady, we assume equilibrium between dissipation and viscous diffusion in the budget equation of T 2 . Combined with the continuity of temperature and heat flux at the fluid-solid interface, one gets:

T 2 f (y) = T 2 f (y f ) + a f (y -y f ) + P r (y -y f ) 2 ε θ,f (y f ) T 2 s (y) = T 2 s (y s ) + a s (y -y s ) + GP r (y -y s ) 2 ε θ,s (y s ) T 2 f (y f s ) = T 2 s (y f s ) a f + 2P r (y f s -y f ) ε θ,f (y f ) = G 2 (a s + 2GP r (y f s -y s ) ε θ,s (y s )) (12)
Once cos (φ) at the fluid-solid interface is obtained -as already explained, we use the value at the first solid cell -it is straightforward to estimate ∂ y T ∂ y T on both sides of the fluid-solid interface. For instance, at the fluid-solid interface and on the fluid side, one gets

∂ y T f ∂ y T f =     ∂ y T f 2 /2 T f 2 cos (φ)     2 (13)
Once the quantities ∂ i T ∂ i T have been extrapolated on both sides of the fluid-solid interface, it is straightforward to estimate the dissipation rates. As the extrapolation procedure proposed here relies deeply on the continuity of the temperature and heat flux, equation ( 9) is automatically satisfied by the reconstructed quantities.

New LES Results

3.1. Case Re τ = 395 and P r = 0.71 In this subsection, the main LES performed is presented, alongside with some key results and a correlation for the discontinuity of ε θ at the fluid-solid interface. As indicated Table 1, the simulation contains 51 passive scalars, 49 being coupled at the fluid-solid interface while the remaining two have either an imposed value or an imposed flux at the fluid boundary. Regarding the 49 coupled scalars, the fluid-solid thermal properties ratio investigated are described in the Table 2. For both the fluid-to-solid thermal diffusivity ratio G i and the thermal activity ratio K j , 7 values are simultaneously investigated. They vary over a range of two decades, centered around unity. The coupled passive scalar with fluid-solid thermal properties ratios G i and K j is denoted as CHT ij .

The left frame of Figure 2 illustrates the impact of the thermal activity ratio K on ε θ . Lower values of K correspond to conjugate cases closer to the isoT one. Oppositely, higher values correspond to conjugate cases closer to the isoQ one. The right frame of Figure 2 shows that statistics, in the fluid domain, are mostly driven by the thermal activity ratio, the impact of the thermal diffusivity ratio G being much weaker. However, inside the solid, the situation is more complex and the impact of G on ε θ is more visible. Overall, the present LES results corroborate well the DNS results obtained at a lower Reynolds number (Re τ = 150) previously published in [START_REF] Flageul | On the discontinuity of the dissipation rate associated with the temperature variance at the fluid-solid interface for cases with conjugate heat transfer[END_REF]. The strong impact of K on the turbulent statistics probably comes from the nature of the flow, which is dominated by wall-normal fluctuations in a turbulent channel, and may not hold for more complex flow configurations. The reconstructed statistics at the fluid-solid interface have a realistic behaviour: as expected, ε θ monotonically decays in the first solid cell.

Proposed correlation

On top of those qualitative observations, the large number of combination of fluid-solid thermal properties ratios investigated allow a more quantitative analysis. We look for power-law correlations for statistics at the fluid-solid interface as follow:

∂ y T f ∂ y T f ∂ i T f ∂ i T f (CHT ij ) = 1 1 + C ε • G xε i • K yε j with C ε > 0 (14)
The main property of equation ( 14) is that the relative wall-normal contribution remains bounded in [0, 1]. The resulting ratio of dissipation rates at the fluid-solid interface will thus remain bounded between K 2 and 1 G , the theoretical bounds. Here, the parameters C ε is derived from the case CHT 44 , which has unit fluid-solid thermal properties ratios. Regarding the coefficients x ε and y ε , they are simply determined with a linear Relative error for the correlation ( 16), in %. regression using logarithms:

x ε log (G i ) + y ε log (K j ) = log ∂ i T f ∂ i T f ∂ y T f ∂ y T f (CHT ij ) -1 -log (C ε ) (15)
In the present study, the linear regression is performed with the Analysis ToolPak available in Microsoft Excel. The resulting correlation at the fluid-solid interface is

ε θ,s ε θ,f = 1 G + K 2 - 1 G 1 1 + 0.0799 • G 0.225 • K 1.90 (16)
The right frame of Figure 3 shows that the correlation for the ratio of ε θ is fairly for the range of thermal properties investigated here: the relative error is below 4 % for most of the cases considered. However, the correlation tends to overestimate the ratio ε θ,s ε θ,f when both K and G are low (lower-left corner), and high (upper-right corner).

Indeed, the analysis proposed here remains qualitative as a better strategy to derive correlations would have been to randomly select 49 combinations of fluid-solid thermal properties ratios inside the interval [0.1, 10]. It seems important to stress that the objective of the present work is not to derive highly accurate correlations for turbulent channel flows as it would probably require a more complex formulation than the one in equation [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF]. The objective is rather to propose a methodology, which can be adapted and applied on any configuration of interest.

Cases at higher Re τ and P r

Regarding the simulations at a higher Re τ and P r, only 3 coupled passive scalars were transported. The left frame of Figure 4 shows the impact of the Reynolds and Prandtl numbers on ε θ for the case with unit ratio of fluid-solid thermal properties. The Reynolds number has a very limited impact on the profile around the fluid-solid interface while the Prandtl number has a more pronounced impact, especially in the fluid domain.

The right frame of Figure 4 shows this impact on ε θ for cases with G = 0.1 and K ≈ 0.23. The case Re τ = 395 and P r = 0.71 has K = 0.2 while the others have K = 0.23. This is representative of pressurized water and steel, except for the value of G, which should be a decade lower. We limited it to 0.1 to remain within the bounds of validity of the correlation.

Overall, the impact of the Reynolds and Prandtl numbers is similar to what was observed at G = K = 1. Remarkably, the case Re τ = 395 and P r = 0.71 exhibits a lower ε θ in the solid domain. On the one hand, this might be the effect of the small difference on K. On the other, it might be a distorsion due to the logarithmic ordinate. More quantitatively, we can measure the relative error produced by the correlation [START_REF] Flageul | Impact of unresolved smaller scales on the scalar dissipation rate in direct numerical simulations of wall bounded flows[END_REF]. The cases with G = K = 1 can be discarded as they all lead to a continuity of ε θ at the fluidsolid interface. The remaining cases are G = 1.3 and K = 2.8 -representative of pressurized water and plexiglas -and G = 0.1 and K = 0.23, representative of pressurized water and steel (except for G, as already explained). For the LES at Re τ = 395 and P r = 0.71, we measure the error of the correlation at the approximate values of (G, K) = (1, 2) and (G, K) = (0.1, 0.2), respectively.

At this stage, it is important to stress that the correlation [START_REF] Flageul | Impact of unresolved smaller scales on the scalar dissipation rate in direct numerical simulations of wall bounded flows[END_REF] has no dependence on Re τ and P r and is thus not expected to perform well at higher Re τ and P r. The amount of data available so far (2 values for Re τ and 2 values for P r) is clearly insufficient to introduce such a dependence. As a workaround, we derive the following correlation from ( 16):

ε θ,s ε θ,f = 1 G + K 2 - 1 G 1 1 + C ε • G 0.225 • K 1.90 (17) with 1 1 + C ε = ∂ y T f ∂ y T f ∂ i T f ∂ i T f (Re τ , P r, G = 1, K = 1)
For the case Re τ = 395 and P r = 0.71, this is exactly the correlation [START_REF] Flageul | Impact of unresolved smaller scales on the scalar dissipation rate in direct numerical simulations of wall bounded flows[END_REF]. For the other cases, it uses the anisotropy of the fluctuating temperature gradient obtained at G = K = 1 to improve the correlation, while keeping the exponents on G and K as in [START_REF] Flageul | Impact of unresolved smaller scales on the scalar dissipation rate in direct numerical simulations of wall bounded flows[END_REF].

As shown Table 3, the modified correlation performs much better compared to the original one. Indeed, this could have been expected for thermal properties ratios close to unity, as in the case (G, K) ≈ (1.3, 2.8). It seems to hold quite well for ratio of thermal properties further away from unity, as in the case (G, K) ≈ (0.1, 0.23).

Discussion and Conclusions

In the present paper, we have developed a strategy to extrapolate statistics at the fluid-solid interface for CFD codes based on finite volume solvers which can be applied to fields coupled across this interface. As it is, the strategy is applicable to wall-resolved LES of statistically steady flows. Extension to wall-modelled LES is not straightforward as the current strategy assumes the first fluid cell is located inside the viscous sublayer.

A correlation for the discontinuity of ε θ at the fluid-solid interface for turbulent channel flows at Re τ = 395 and P r = 0.71 is also proposed, as expressed in equation ( 16). The correlation is valid for fluid-solid thermal properties ratios G and K within [0.1, 10]. Extension to thermal properties ratios further away from unity could be achieved but special care should be taken. First, transient might become excessively long and statistical convergence in the solid domain might be hard to reach when G becomes too large. In addition, as the coupled case degenerates towards an imposed temperature or an imposed heat flux, temperature fluctuations in the solid domain might get damped over a very short distance and a very short time scale, thus requiring a reduced time step and a very fine mesh in the solid domain.

A modification of the correlation towards cases at a higher Reynolds number and different Prandtl number, as expressed in equation [START_REF] Kasagi | Direct numerical simulation of passive scalar field in a turbulent channel flow[END_REF], is proposed and assessed with three LES. The strategy is promising but would require more simulations to be thoroughly validated. Regarding conjugate heat transfer in turbulent flows, the anisotropy of the fluctuating temperature gradient at the fluid-solid interface clearly remains the key quantity, both for the proposed correlations and for future models. This work suggests that, given a flow regime, the impact of the fluidsolid thermal properties ratios can be clearly isolated. Indeed, much more work remains to be done to confirm this point.

kl 5/ 4 (

 4 S kl S kl ) 5/2 + S d kl S

Figure 1 :

 1 Figure 1: Channel flow with fluid-solid thermal coupling. The solid domains are on top (y > δ) and bottom (y < -δ) of the fluid one (-δ < y < δ). Dotted lines: fluid-solid interfaces at y = ±δ.

Figure 2 :

 2 Figure 2: Scalar dissipation rate for the wall-resolved LES at Re τ = 395 and P r = 0.71. Left frame: Cases G = 0.2. Right frame: Cases K = 0.2, with a focus on the near-wall behaviour in the fluid. y + < 0: Solid domain. y + > 0: Fluid domain.

Figure 3 :

 3 Figure 3: Contour plots at the fluid-solid interface for the wall-resolved LES at Re τ = 395 and P r = 0.71. Left frame: Ratio of the solid and fluid scalar dissipation rates. Right frame: Relative error for the correlation (16), in %.

Figure 4 :

 4 Figure 4: Scalar dissipation rate for the wall-resolved LES at various Re τ and P r. Left frame: Cases G = 1 and K = 1. Right frame: Cases G = 0.1 and K ≈ 0.23.

Table 1 :

 1 LES performed

	L x /δ	6.283	6.283	L x /δ
	L z /δ	3.142	3.142	L z /δ
	δ S /δ	0.375	0.147	δ S /δ
	Re τ	395	1020	Re τ
	P r 0.71 1 0.71 1 P r
	Scalars 51 5	5	5 Scalars

Table 2 :

 2 Thermal properties ratios investigated in the subsection 3.1. . Solid-to-fluid thermal conductivity ratio G 2 .

				K j
		0.1 0.2 0.5	1	2	5	10
	0.1 32 16 6.3 3.2 1.6 0.63 0.32
	0.2 22 11 4.5 2.2 1.1 0.45 0.22
	0.5 14 7.1 2.8 1.4 0.71 0.28 0.14
	G i 1	10 5	2	1	0.5	0.2	0.1
	2 7.1 3.5 1.4 0.71 0.35 0.14 0.071
	5 4.5 2.2 0.89 0.45 0.22 0.089 0.045
	10 3.2 1.6 0.63 0.32 0.16 0.063 0.032

Table 3 :

 3 Relative error in % for the correlations (16) and (17) Correlation (16) (G, K) ≈ (1.3, 2.8) 2.19 6.49 8.69 12.7 Correlation (17) (G, K) ≈ (1.3, 2.8) 2.19 1.21 1.71 1.01 Correlation (16) (G, K) ≈ (0.1, 0.23) 11.7 13.8 15.4 17.7 Correlation (17) (G, K) ≈ (0.1, 0.23) 11.7 8.02 7.78 5.00

	Re τ	395		1020
	P r 0.71	1	0.71	1

https://github.com/code-saturne/code_saturne

Acknowledgements

The authors thank the Slovenian Research Agency and EDF R&D for funding the study (research projects P2-0026 and PR-07184).