C Flageul
email: cedric.flageul@ijs.si

I Tiselj
email: iztok.tiselj@ijs.si

S Benhamadouche
email: sofiane.benhamadouche@edf.fr

M Ferrand
email: martin.ferrand@edf.fr

QUANTIFICATION OF THE DISCONTINUITY OF THE TEMPERATURE VARIANCE DISSIPATION RATE: WALL-RESOLVED LES OF TURBULENT CHANNEL FLOW WITH CONJUGATE HEAT TRANSFER

Keywords: Conjugate heat transfer, wall-resolved LES, Temperature variance dissipation rate

). Actually, there is currently no coupled RANS model for conjugate heat transfer taking this discontinuity into account. As a result, from an industrial perspective, LES remains the best option for thermal fatigue prediction but needs refinement at the wall, which makes it very expensive if not unaffordable, at high Reynolds numbers. In this paper, we will assess the ability of wall-resolved LES to estimate this discontinuity of ε θ on channel flows ( ) using Code_Saturne, EDF in-house and open-source CFD software.

INTRODUCTION

Conjugate heat transfer describes the thermal coupling between a fluid and a solid. It is of prime importance in industrial applications where fluctuating thermal stresses are a concern. As far as the nuclear industry is concerned, relevant applications include severe emergency cooling (Pressurized Thermal Shock) or long-term ageing of materials (T junctions). For such complex applications, characterized by high Reynolds number, investigations often rely on experiments, high Reynolds RANS or wall-modelled LES. However, experimental data on conjugate heat transfer are scarce as walls in lab rigs are often made of plexiglas and the transported scalar studied is often a dye. The development of RANS models for conjugate heat transfer is relatively recent (Craft et al., [START_REF] Craft | Towards the development of RANS models for conjugate heat transfer[END_REF]).

In [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][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 authors have established that in case of conjugate heat-transfer, the dissipation rate (ε θ ) associated with the halved variance of the temperature (θ' 2 /2) is discontinuous at the fluid-solid interface; thus, challenging the ability of (U)RANS models to investigate configurations with conjugate heattransfer. As a result, from an industrial perspective, LES remains the best option for thermal fatigue prediction but needs refinement at the wall, which makes it very expensive if not unaffordable, at high Reynolds numbers.

Experimental data on conjugate heat transfer are scarce, expensive and difficult to carry out in order to obtain fine quantities such as ε θ . This makes DNS and wall-resolved LES the only solutions to obtain such data. In this paper, in order to deal with higher Reynolds numbers than in DNS, we will assess the ability of wall-resolved Large Eddy Simulation to estimate this discontinuity of ε θ on channel flows using Code_Saturne, EDF in-house and open-source CFD software. This is a step forward towards a rich validation database for future RANS models adapted to conjugate heat transfer.

The structure of this paper is as follows. In the second section, governing equations are presented alongside with their numerical discretization. In the third section, the cases simulated are fully described. In the fourth section, the results are presented. In the fifth section, they are further discussed alongside with their consequences for RANS and LES modeling.

GOVERNING EQUATIONS AND DISCRETIZATION

The equations presented hereafter assume that the fluid is Newtonian with constant physical properties and that the flow is incompressible. They also assume that the thermal properties in the solid domain are constant, although they may differ from the fluid ones. Naturally, the discontinuity of ε θ at the fluid-solid interface holds for compressible flows or when physical properties are variable.

Governing Equations

The continuity equation reads:

(1)

The momentum equation for direction reads: with

The source term is a constant pressure gradient and is active only in the streamwise direction. When using the WALE model, the subgrid-scale (SGS) turbulent viscosity is defined by: with

In equation ( 3), is the traceless symmetric part of the square of the velocity gradient tensor, as defined by Nicoud and Ducros, [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF]. In case of conjugate heat transfer, the energy equation reads:

in the fluid domain . ( 4
)
in the solid domain .

(

) 5
and at the fluid-solid interface .

In equation ( 4), the turbulent Prandtl number , is set to (Grötzbach [START_REF] Grötzbach | Revisiting the resolution requirements for turbulence simulations in nuclear heat transfer[END_REF]) and the source term is as defined by Kasagi et al., [START_REF] Kasagi | Direct numerical simulation of passive scalar field in a turbulent channel flow[END_REF]. In equation [START_REF] Kasagi | Direct numerical simulation of passive scalar field in a turbulent channel flow[END_REF], is the derivative along the vector normal to the fluid-solid interface. This equation is a direct transcription of the continuity of the temperature and heat-flux across the interface. Following Tiselj and Cizelj, [START_REF] Tiselj | DNS of turbulent channel flow with conjugate heat transfer at Prandtl number 0.01[END_REF], 2 independent dimensionless numbers describe this fluid-solid thermal coupling: and

In equation [START_REF] Tiselj | DNS of turbulent channel flow with conjugate heat transfer at Prandtl number 0.01[END_REF], is the fluid-to-solid thermal diffusivity ratio and is the thermal activity ratio, which is also the fluid-to-solid thermal effusivity ratio. The thermal properties ratios (alongside with the Reynolds and Prandtl number) used in our simulations are further described in section 3.

Following Flageul et al., [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][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], in case of conjugate heat transfer, ε θ is discontinuous at the fluid-solid interface and satisfies:

at the fluid-solid interface .

As equation ( 8) is a convex combination, the ratio of the solid-to-fluid dissipation rates at the fluid-solid interface is bounded by and . Equation ( 8) holds as long as the vector normal to the fluid-solid interface is well-defined, which means any flat or curved interface. It becomes ill-defined only at edges or corners.

In [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][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], equation [START_REF] Archambeau | Code_Saturne: a finite volume code for the computation of turbulent incompressible flows -industrial applications[END_REF] was derived in the absence of SGS model ( ) using dissipation rates defined by: and

In this paper, we will examine the validity of equation ( 8) for wall-resolved LES of the turbulent channel flow using the WALE model, a constant turbulent Prandtl number of 0.5 and dissipation rates defined by: and

, the dissipation rate of the turbulent kinetic energy is similarly defined by: (11)

Discretization

Our LES are performed with the open source solver Code_Saturne (www.code-saturne.org), an unstructured collocated finite volume solver for incompressible flows. It uses a SIMPLEC algorithm for pressure-velocity coupling, which requires Rhie and Chow interpolation to avoid odd-even oscillations.

Since version 5.0, Code_Saturne can be used stand-alone to perform conjugate heat-transfer with simultaneous solving of the fluid and solid temperature.

As our meshes are made of orthogonal hexahedra, no gradient reconstruction sweeps are needed. In addition, the existing slope-test by default is deactivated for the scalars. As a result, for the velocity and the temperature, a fully centered scheme is used, combined with a Crank-Nicolson time scheme. The convective term is linearized and the transporting velocity uses an Adams-Bashforth time scheme. Further details about Code_Saturne 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].

CASE SETUP

The structure of this section is as follows. In the first subsection, the validation case is presented. It corresponds to a Reynolds number of based on the bulk velocity ( based on the friction velocity) and uses DNS data from Flageul et al. [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][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]. In the second subsection, 2 other cases are described. They both correspond to a Reynolds number of based on the bulk velocity ( based on the friction velocity) but the thickness of the solid part differs.

In the present study, x, y and z stand for the streamwise, wall-normal and spanwise directions, respectively. Periodic boundary conditions are used in the directions x and z. For all the configurations, the extension of the solid domain in the directions x and z is identical to the fluid one.

All the DNS and LES are performed with a Prandtl number ( ) of 0.71. The Reynolds number ( ) based on the bulk velocity is defined with the kinematic viscosity , the channel half-height and the averaged bulk velocity, the two latter equaling unity in the present study. For cases with conjugate heat transfer, an imposed heat flux is always used at the outer wall.

For each case, 2 LES are performed: one without SGS model ( ) and one with the WALE model combined with a constant turbulent Prandtl number ( ) of 0.5.

Validation Case at Re=2280

For the validation case at Re=2280, Table I summarizes the simulation parameters. 11 passive scalars are transported by the flow. One passive scalar satisfies an imposed temperature at the fluid-wall (isoT), another satisfies an imposed heat flux at the fluid-wall (isoQ) and the 9 remaining passive scalars are subjected to fluid-solid coupling with thermal properties ratios around unity, as in Flageul et al. [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][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].

Cases at Re=7060

For the cases at Re=7060, Table II summarizes the simulation parameters. For the first case (case 1), the solid domain is as thick as in the validation case and we have performed a DNS with the pseudo-spectral code previously used by Tiselj and Cizelj, [START_REF] Tiselj | DNS of turbulent channel flow with conjugate heat transfer at Prandtl number 0.01[END_REF]. Case 1 and the DNS contain only 3 scalars. The first satisfies an imposed temperature at the fluid-wall (isoT), the second an imposed heat flux at the fluid-wall (isoQ) and the third fluid-solid coupling with the same thermal properties in the fluid and in the solid.

For the second case (case 2), the solid domain is thinner (values between brackets in table II), thus allowing us to reduce the duration of our simulation. In wall-units, the solid domain remains thick enough so that there is no impact of the boundary condition used at the outer wall on the statistics at the fluidsolid interface. As the duration of the simulation is reduced, we can transport more scalars while keeping the computational cost of the simulations within reasonable bounds.

For case 2, 51 scalars are transported. One satisfies an imposed temperature at the fluid-wall (isoT), the second an imposed heat flux at the fluid-wall (isoQ) and the 49 others fluid-solid coupling. 7 values for the thermal diffusivity ratio and the thermal activity ratio are examined. For both parameters and , the same range was examined, spanning 2 decades and centered around unity: .

RESULTS

Validation Results at Re=2280

LES models are designed for flows with a relatively high Reynolds number. As this validation case uses a relatively low Reynolds number, only qualitative agreement is expected.

Averaged velocity, Reynolds stresses and dissipation rate ε at Re=2280

The averaged velocity and dissipation rate ε in Figure 1 show a relatively good agreement with DNS. The agreement is more qualitative for the Reynolds stresses, as shown in Figure 2.

Scalar-related statistics at Re=2280

As the agreement with DNS is similar for most of the scalars, only 3 are presented here: isoT, isoQ, and [G,K]=[0.5, ]. Hereafter, the fluid (solid) domain is located at y + >0 (y + <0). In figure 3, the temperature variance shows a qualitative agreement with the DNS. The wall-normal contribution in ε θ , which drives the discontinuity (following equation ( 8)) is slightly overestimated at the fluid-solid interface for the conjugate cases. In figure 4, the discontinuity of ε θ at the fluid-solid interface is well estimated by the present LES, although the overall agreement in the fluid domain is only qualitative.

Case 1 at Re=7060

Averaged velocity, Reynolds stresses and dissipation rate ε, case 1 at Re=7060

The averaged velocity and dissipation rate ε in Figure 5 show a relatively good agreement with DNS. The agreement is still qualitative for the Reynolds stresses, as shown in Figure 6. Overall and as expected, LES performed at Re=7060 are closer to DNS, compared with the validation case at Re=2280.

Scalar-related statistics, case 1 at Re=7060

In figure 7, the temperature variance shows a qualitative agreement with the DNS. Our LES compare favorably with the DNS in the solid domain and at the fluid-solid interface. For the dissipation rates ε θ plotted Figure 8, our LES results show reasonable agreement with the DNS data, especially at the fluidsolid interface where the WALE model outperforms the simulation without SGS model. For the relative wall-normal contribution in ε θ , Figure 9 suggests that our LES overestimate it at the fluid-solid interface for conjugate cases, as observed at a lower Reynolds number (Figure 3).

Case 2 at Re=7060

For case 2, average velocity, Reynolds stresses and dissipation rate ε are not plotted as they are similar to the ones in case 1.

Scalar-related statistics, case 2 at Re=7060

Due to the very large number of scalars transported, only partial results supporting the main conclusions are presented here. Figure 10 shows that in the fluid domain, the variance of the temperature depends mostly on the thermal activity ratio K. It is similar for the dissipation rate, which has a more complex behavior in the solid domain (Figure 11), which clearly shows the discontinuity of ε θ . This discontinuity is driven by the relative wall-normal contribution in ε θ , plotted in Figure 12.

DISCUSSION

Qualitatively, the results show that our wall-resolved LES can quantify the discontinuity of ε θ at the fluidsolid interface. Using equations ( 8) and (10), the error in our statistics can be evaluated. However, our statistics are defined at the cell center. As a result, we have no direct access to the quantities present in equation ( 8) as the fluid-solid interface is located at faces of cells. Fortunately, our LES are wall-resolved, so those quantities can be reconstructed with a relatively small error (first order in space). First, we recall equation [START_REF] Archambeau | Code_Saturne: a finite volume code for the computation of turbulent incompressible flows -industrial applications[END_REF] and reorder it to use the relative wall-normal contribution in the solid domain:

at . ( 12
) at . ( 13
)
However, following figure 12, one may say that the relative wall-normal contribution at y + =±0.5 is only a crude approximation of the relative wall-normal contribution at the fluid-solid interface. Using Taylor expansion, the variance of the temperature at the fluid-solid interface can be reconstructed more accurately. We reconstruct it with a quadratic Taylor expansion. We also reconstruct the relative wallnormal contribution at the fluid-solid interface using a linear Taylor expansion. This is further described in appendix A. Table III gives the value of the square root of the product of equations ( 12) and ( 13) for the WALE simulation. Results are similar for the simulation without SGS model. Values written in a bold font and underlined indicate an error above 10%. For most of the scalars considered here, the combination of the statistical error and the reconstruction approximation is limited to a few percent, or less. However, cases with a high K and a high G give rise to higher errors. The source of this higher error is currently unknown. High values of G are associated with large time scales for the temperature fluctuations in the solid domain. Therefore, the present results at high G and K may be biased by statistical convergence and need deeper analysis.

At the fluid-solid interface, θ' 2 is less than 10% of the isoQ value for all the conjugate cases with K≤0.2, and for the cases with K=0.5 and G≤0.2. From this perspective, the present conjugate cases quickly converge towards the isoT case for thermal activity ratios below unity. The maximal value of θ' 2 is around 84% of the isoQ one and is obtained with K=10 and G=10. One may conclude that the thermal properties ratios studied here are not far enough from unity to reach the isoQ case. However, as already stated, the present results at high G and K may be biased by statistical convergence and need deeper analysis.

The present simulations resulted in a relatively huge dataset. The authors are currently analyzing it. Our main objective is to extract scaling laws from the statistics related to the coupled scalars at the fluid-solid interface and in the solid domain, with the hope that those scaling laws will be useful for future (U)RANS turbulence models and LES wall-models.

CONCLUSIONS

The results presented here show that wall-resolved LES can be used to quantify the discontinuity of ε θ at the fluid-solid interface in case of conjugate heat transfer. The combination of the WALE SGS model and a constant turbulent Prandtl number of 0.5 showed good agreement with available DNS results. Using equations ( 12) and ( 13), we also propose an indicator allowing one to evaluate the error at the fluid-solid interface. This indicator is useful as it is stand-alone and derived solely from the continuity of temperature and heat flux at the fluid-solid interface.

From the perspective of a rich validation database for future (U)RANS models for conjugate heat transfer adapted to the discontinuity of ε θ , this is encouraging as more complex configurations, closer to industrial applications, can be investigated on the basis of the present work, which uses a relatively simple SGS model. The proposed error indicator will give an objective criterion to assess the accuracy of the statistics in the aforementioned future database.

However, much more work remains necessary to confidently assess industrial challenges. Aside from more complex geometries, an interesting perspective is the extension towards active scalar cases for which the fluid-solid coupling can potentially change the whole flow and heat transfer. This will require a careful examination to make sure the present SGS model handles correctly buoyancy effects or variable fluid properties, such as the ones encountered in supercritical applications.

Figure 1 .

 1 Figure 1. Left: averaged streamwise velocity. Right: estimate of the total turbulent dissipation ε. Black solid line: DNS. Red dashed line: No SGS. Blue dotted line: WALE.

Figure 2 .

 2 Figure 2. Reynolds stresses. Black solid line: DNS. Red dashed line: No SGS. Blue dotted line: WALE.

Figure 3 .

 3 Figure 3. Left: temperature variance. Right: wall-normal contribution in ε θ,f . Black: DNS. Red: No SGS. Blue: WALE. Dashed lines: isoT. Dotted lines: isoQ. Solid lines: G=0.5, K= .

Figure 4 .

 4 Figure 4. Left and right: ε θ . Black: DNS. Red: No SGS. Blue: WALE. Dashed lines: isoT. Dotted lines: isoQ. Solid lines: G=0.5, K= .

Figure 5 .Figure 6 .

 56 Figure 5. Left: averaged velocity. Right: ε. Black solid line: DNS. Red dashed line: No SGS. Blue dotted line: WALE.

Figure 7 .

 7 Figure 7. Left and right: temperature variance. Black: DNS. Red: No SGS. Blue: WALE. Dashed lines: isoT. Dotted lines: isoQ. Solid lines: G=1, K= .

Figure 8 .

 8 Figure 8. Left and right: ε θ . Black: DNS. Red: No SGS. Blue: WALE. Dashed lines: isoT. Dotted lines: isoQ. Solid lines: G=1, K= .

Figure 9 .

 9 Figure 9. Left and right: wall-normal contribution in ε θ,f . Black: DNS. Red: No SGS. Blue: WALE. Dashed lines: isoT. Dotted lines: isoQ. Solid lines: G=1, K= .

Figure 10 .

 10 Figure 10. Left and right: temperature variance. Left: G=0.2. Right: K=0.2. Black: DNS. Red: No SGS. Blue: WALE. Dashed lines: isoT. Dotted lines: isoQ. Solid lines: Conjugate cases.

Figure 11 .

 11 Figure 11. ε θ . Left: G=0.2. Right: K=0.2. Black: DNS. Red: No SGS. Blue: WALE. Dashed lines: isoT. Dotted lines: isoQ. Solid lines: Conjugate cases.

Figure 12 .

 12 Figure 12. Wall-normal contribution in ε θ,f . Left: G=0.2. Right: K=0.2. Black: DNS. Red: No SGS. Blue: WALE. Dashed lines: isoT. Dotted lines: isoQ. Solid lines: Conjugate cases.

Table I . Simulation parameters for the case Re=2280

 I

		Present	Flageul et al. [2-3]
	Fluid domain extension	[12.56, 2, 6.28]	[25.6, 2, 8.52]
	Solid thickness	1	1
	Grid	[64, 62, 64]	[256, 193, 256]
		[1, 15]	[0.49, 4.8]
		[30, 15]	[14.8, 5.1]
	Time step (wall-units)	0.25	0.02
	Duration (wall-units)	25'000	29'000

Both the reference DNS and the present LES use long averaging time to ensure statistical convergence deep inside the solid domain. The thick solid domain used prevents the boundary condition used at the outer wall from interfering with the statistics at the fluid-solid interface.

Table II . Simulation parameters for the cases at Re=7060

 II

		Present	DNS
	Fluid domain extension	[6.283, 2, 3.142]	[6.283, 2, 3.142]
	Solid thickness	1 (0.375)	1
	Grid	[81, 106, 81]	[256, 257, 256]
		[1, 15]	[0.03, 4.8]
		[30, 15]	[9.6, 4.8]
	Time step (wall-units)	0.22	0.016
	Duration (wall-units)	220'000 (132'000)	9'500

Table III . Geometrical mean of the reconstructed dissipation ratios in equations (12) and (13), case 2, WALE simulation. Theoretical value is 1.

 III

		K=0.1	K=0.2	K=0.5	K=1	K=2	K=5	K=10
	G=0.1	1.002	1.005	1.011	1.018	1.022	1.015	1.000
	G=0.2	1.001	1.003	1.007	1.010	1.010	1.000	0.984
	G=0.5	1.000	1.000	1.001	1.001	1.000	0.995	0.991
	G=1	0.999	0.999	0.999	1.000	1.002	1.009	1.022
	G=2	0.999	0.999	1.000	1.004	1.012	1.039	1.079
	G=5	0.999	1.000	1.005	1.015	1.039	1.108	1.206
	G=10	1.000	1.002	1.011	1.030	1.073	1.193	1.363

ACKNOWLEDGMENTS

The authors thank EDF R&D and the institute Jožef Stefan for funding the present study and providing computational time.

APPENDIX A

In the first cell, both in the fluid and in the solid domain, we approximate the variance of the temperature with a quadratic Taylor expansion: in the first fluid cell.

(14a) in the first solid cell.

(14b)

The center of the first cells are located at and while the fluid-solid interface is located at . The variance of the temperature at the cell centers ( and ) is known. We assume the flow is statistically steady. In addition, in the first cell, we assume equilibrium of dissipation and molecular diffusion in the budget equation of , which is reasonable for a wall-resolved LES. This leads to and (15)

As the dissipation rates at the cell centers are known, 2 unknown remains. They are derived from the continuity of temperature and heat flux:

The ratio of thermal conductivities can be derived from the thermal properties ratios:

The resulting 2x2 system is:

.

We have injected the solution in equation (14a) to compute .

The terms and are treated similarly using a linear Taylor expansion, their continuity at the fluid-solid interface and the discontinuity of their wall-normal derivative, which behaves like the derivative of the temperature variance at the fluid-solid interface.

At the fluid-solid interface, the term is also discontinuous:

However, for this last quantity, no compatibility condition involving its derivative is available so we don't use a linear Taylor expansion. We use the normalized one-point correlation associated with the temperature and its wall-normal derivative, which is continuous at the fluid-solid interface:

at the fluid-solid interface.

(20)

On both sides of equation ( 20), the numerator is the halved derivative of the temperature variance. This quantity and the root mean square of the temperature are available at the fluid-solid interface using equations (14a) and (14b). In addition, we assume that this normalized correlation at the fluid-solid interface is the arithmetic mean of the corresponding values at the first fluid and solid cell center. This strategy allows us to reconstruct the relative wall-normal contribution at the fluid-solid interface.


Cédric Flageul

Iztok Tiselj

Impact of unresolved smaller scales on the scalar dissipation rate in direct numerical simulations of wall bounded flows

Keywords: Passive scalar, Dissipation rate, Anisotropy, Direct Numerical Simulation, Channel flow, Sampling error

. All scalar fields are computed in parallel using a single velocity field resolved with the finest resolution, thus reducing the statistical variability.

In addition, to confidently assess the grid spacing requirement, we also evaluate the statistical uncertainty. The "standard" resolution of the DNS studies (resolution used by Kim et al. [2]) is usually sufficient for predictions of first and second-order integral turbulence scalar field statistics. Non-negligible corrections of the fourth-order integral statistics, especially the scalar dissipation variance profile, are observed with enhancement of the scalar resolution from the one used in the "standard" DNS studies to the resolution recommended by Vreman and Kuerten [3], which is roughly two times finer in each spatial direction. Further resolution enhancements produce only marginal differences.

Introduction

Since the pioneering efforts of Obukhov [START_REF] Obukhov | Structure of the temperature field in turbulent flow[END_REF], Corrsin [START_REF] Corrsin | On the spectrum of isotropic temperature fluctuations in an isotropic turbulence[END_REF] and Batchelor et al.

( [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid part 1. general discussion and the case of small conductivity[END_REF][START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid part 2. the case of large conductivity[END_REF]), passive scalars in turbulent flows have been the focus of a number of studies. As summarized in the review by Warhaft [START_REF] Warhaft | Passive scalars in turbulent flows[END_REF], experiments and simulations are challenging classical descriptions of passive scalars derived from the Kolmogorov cascade phenomenology. Evidences suggest a strong coupling between large and small scales and no local isotropy at inertial and dissipation scales. They also suggest that passive scalars are associated with a stronger intermittency compared with the velocity.

Following Batchelor et al. ([6,[START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid part 2. the case of large conductivity[END_REF]), for a unit Prandtl number as herein, the smallest spatial scale for scalar mixing is the Kolmogorov scale η = ν 3 /ε 1/4 , where ν is the dynamic viscosity and ε the mean dissipation rate of the turbulent kinetic energy. Similarly, the smallest time scale is the Kolmogorov time scale τ η = ν/ε. As the passive scalar is intermittent, locally, structures with a spatial (temporal) span shorter than η (τ η ) appear in a flow. Suspecting those fine structures to be related to highly dissipative events, a number of DNS have been performed with sub-Kolmogorov scales resolved (Schumacher et al. [START_REF] Schumacher | Very fine structures in scalar mixing[END_REF], Donzis and Yeung [START_REF] Donzis | Resolution effects and scaling in numerical simulations of passive scalar mixing in turbulence[END_REF], Galantucci and Quadrio [START_REF] Galantucci | Very fine near-wall structures in turbulent scalar mixing[END_REF]).

One of the central quantity in those studies is ε θ , the scalar dissipation rate, defined by

ε θ = 2α|∇θ| 2 = 2α 3 i=1 (∂ i θ) 2 (1)
where α is the thermal diffusivity. According to Pope [START_REF] Pope | Small scales, many species and the manifold challenges of turbulent combustion[END_REF], this quantity he calls the all-important dissipation rate matters in combustion models. It is also a central quantity in Reynolds Averaged Navier-Stokes (RANS) turbulence models as it appears in the budget equation of the scalar variance. Lately, 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] showed that the dissipation rate associated with the temperature variance is discontinuous at the fluid-solid interface in case of conjugate heat transfer. This is prominent for industrial applications where thermal fatigue is a concern.

Studying homogeneous isotropic turbulence, Schumacher et al. [START_REF] Schumacher | Very fine structures in scalar mixing[END_REF] showed that the improved resolution matters when investigating the tails of the Probability Density Function (PDF) of ε θ , which correspond to low probability events associated with high or low dissipation rates. More specifically, they show that a poorer resolution has a stronger impact on regions of low ε θ than on those of high ε θ . On a similar configuration, Donzis and Yeung [START_REF] Donzis | Resolution effects and scaling in numerical simulations of passive scalar mixing in turbulence[END_REF] showed accurate estimation of advanced statistics (scalar dissipation intermittency exponent, structure functions at moderately high orders and PDF of ε θ up to 200 ε θ ) with a grid spacing equal to the Batchelor scale, which is exactly the Kolmogorov scale in the present study.

Studying a turbulent channel flow (Re τ = 160, P r = 1), Galantucci and

Quadrio [START_REF] Galantucci | Very fine near-wall structures in turbulent scalar mixing[END_REF] extended the analysis to wall-bounded flows. They showed resolution effects on the profiles of the mean ε θ and its variance, but also on the PDF of ε θ and recommended a very fine spatial resolution (∆ + x = ∆ + z = 1 and 0.43 < ∆ + y < 2). In the streamwise (spanwise) direction, this is 6 (4) times finer than what is necessary for the velocity according to Vreman and Kuerten [START_REF] Vreman | Comparison of direct numerical simulation databases of turbulent channel flow at re τ = 180[END_REF]. The authors of the present study estimate that the resolution proposed by Vreman and Kuerten is sufficient for accurate predictions of the key integral turbulent statistics of the passive scalar field at P r = 1, including the average scalar dissipation rate and its variance. The main objective of the present paper is to assess this claim.

The structure of the paper is as follows. In the second section, the governing equations and the computational setup are described alongside with the procedure to estimate the sampling error. In the third section, preliminary investigation on coarser grids is presented. In the fourth section, the DNS results are presented. Discussion and conclusions are collected in the last section.

Governing equations, computational setup and sampling error

Dimensionless equations of the incompressible turbulent channel flow with transport of a passive scalar can be found in various sources (Kasagi et al. [START_REF] Kasagi | Direct numerical simulation of the passive scalar field in a two-dimensional turbulent channel flow[END_REF], Kawamura et al. [START_REF] Kawamura | Dns of turbulent heat transfer in channel flow with low to medium-high prandtl number fluid[END_REF]):

∇.u = 0
(2)

∂ t u = -∇. (u.u) + 1 Re τ ∇ 2 u -∇ (p) + -→ 1 x (3)
∂ t θ = -∇. (u.θ) + 1 Re τ P r ∇ 2 θ (4)
Equations ( 3) and ( 4) are normalized with the channel half width h, the kinematic viscosity ν and the friction velocity u τ . Low friction Reynolds number Re τ = 160 and Prandtl number P r = 1 were selected in order to replicate the conditions of the simulations performed by Galantucci and Quadrio [START_REF] Galantucci | Very fine near-wall structures in turbulent scalar mixing[END_REF]. As the Prandtl number is unity, Kolmogorov and Batchelor length-scales are identical.

Periodic boundary conditions are used in the streamwise and spanwise directions, labelled x and z, respectively, while the wall-normal direction is labelled y. The forcing term -→ 1 x represents a constant pressure gradient in the streamwise direction and has a unit amplitude thanks to the normalization used. Boundary conditions for the passive scalar fields at the channel walls are θ = 1 at y = 1, and θ = -1 at y = -1 and were previously used in the simulations of Papavassiliou and Hanratty [START_REF] Papavassiliou | Transport of a passive scalar in a turbulent channel flow[END_REF], Johansson and Wikström [START_REF] Johansson | Dns and modelling of passive scalar transport in turbulent channel flow with a focus on scalar dissipation rate modelling[END_REF], and Galantucci

and Quadrio [START_REF] Galantucci | Very fine near-wall structures in turbulent scalar mixing[END_REF].

The equations are solved with a pseudo-spectral scheme. Fourier series are used in the x and z directions and Chebyshev polynomials are used in the y direction. Second-order accurate time differencing (Crank-Nicolson scheme for diffusive terms and Adams-Bashforth scheme for other terms) is used with maximum Courant number kept at approximately 0.1. The aliasing error is removed with computation of the nonlinear terms on a grid 1.5 times finer in all directions. The computer code is based on the code developed by Gavrilakis et al. [START_REF] Gavrilakis | Direct and large eddy simulation of turbulence, notes on numerical fluid mechanics[END_REF], which was later modified by Lam and Banerjee [START_REF] Lam | Investigation of turbulent flow bounded by a wall and a free surface[END_REF]. The code was used and verified in simulations of Tiselj et al. ([19], [START_REF] Tiselj | Dns of turbulent channel flow with conjugate heat transfer at prandtl number 0.01[END_REF], [START_REF] Tiselj | Tracking of large-scale structures in turbulent channel with direct numerical simulation of low prandtl number passive scalar[END_REF]).

The extension of the computational domain, normalized with the channel half width h, was taken from the work of Galantucci and Quadrio [START_REF] Galantucci | Very fine near-wall structures in turbulent scalar mixing[END_REF]: L y = 2, and L z = 2.09. Both in the streamwise and spanwise directions, this is about 3 times smaller than the domain used in Kasagi et al. [START_REF] Kasagi | Direct numerical simulation of the passive scalar field in a two-dimensional turbulent channel flow[END_REF]. Such a small computational domain neglects an important part of the large scale structures in the turbulent flow, however, it is known to be sufficient for special studies focused on small scale turbulent structures of the velocity field (Jimenez and Moin [START_REF] Jiménez | The minimal flow unit in near-wall turbulence[END_REF]). The small domain offers a platform for simplified studies of the resolution requirements. It is often overlooked that in addition to obliterating large scale structures, small domain can be affected by a significant sampling error. For instance, Galantucci and Quadrio [START_REF] Galantucci | Very fine near-wall structures in turbulent scalar mixing[END_REF] report up to 5% differences in the friction temperature values in their simulations performed on different resolutions and averaged over the time interval of 2400 viscous time units (statistics based on 60 instantaneous fields).

L x = 4.
In the present study, the resolution requirement for the passive scalar field is examined with a single DNS run. The velocity field is calculated on the finest grid of N x * N y * N z = 360 * 181 * 180 modes (Case 6S in Table 1). Six distinct passive scalar fields are simultaneously transported by this velocity field and resolved with different number of modes, see Table 1. The naming scheme for the passive scalar fields is inspired by the name of the authors who promoted certain resolution. For instance, in 1KMM, KMM stands for Kim, Moser and Moin [START_REF] Kim | Turbulence statistics in fully developed channel flow at low reynolds number[END_REF]. Similarly, in 4VK, VK stands for Vreman and Kuerten [START_REF] Vreman | Comparison of direct numerical simulation databases of turbulent channel flow at re τ = 180[END_REF]. The number preceding the letters corresponds to a disputable ranking of the grids:

1 for the coarsest, 6 for the finest. For the five scalar fields resolved with a lower resolution, all Fourier and Chebyshev modes above the indicated resolution are set to zero at the end of each time step.

This approach, which is comparable to the one in Brethouwer et al. [START_REF] Brethouwer | Micro-structure and lagrangian statistics of the scalar field with a mean gradient in isotropic turbulence[END_REF] or Gotoh et al. [START_REF] Gotoh | Spectral compact difference hybrid computation of passive scalar in isotropic turbulence[END_REF], reduces the statistical variability and eases the separation of the sampling error from the error induced by a coarser spatial resolution.

The separation of statistical uncertainty and resolution effects is of particular importance for the present work: the smaller extension of the domain in the homogeneous directions increases the sampling error, which easily exceeds the tiny differences induced by the variable resolution of the passive scalar field, except for simulations with a very long duration.

The finer resolution (case 6S in Table 1) used in the wall-normal direction follows Vreman and Kuerten [START_REF] Vreman | Comparison of direct numerical simulation databases of turbulent channel flow at re τ = 180[END_REF]. In the streamwise and spanwise directions, case 6S corresponds to the resolution used by Galantucci and Quadrio [START_REF] Galantucci | Very fine near-wall structures in turbulent scalar mixing[END_REF] in their "Medium" simulation. This is 2 to 3 times finer than the recommendation of Vreman and Kuerten [START_REF] Vreman | Comparison of direct numerical simulation databases of turbulent channel flow at re τ = 180[END_REF]. Cases 5GQM and 2GQC correspond to the "Medium"

and "Coarse" cases in Galantucci and Quadrio [START_REF] Galantucci | Very fine near-wall structures in turbulent scalar mixing[END_REF], respectively. The case 4VK is using the resolution recommended by Vreman and Kuerten [START_REF] Vreman | Comparison of direct numerical simulation databases of turbulent channel flow at re τ = 180[END_REF]. The case 3S is similar to 4VK except that it uses a coarser wall-normal grid. Lastly, the resolution in case 1KMM is comparable with most of the previous DNS simulations (Kim et. al. [2], Kasagi et al. [START_REF] Kasagi | Direct numerical simulation of the passive scalar field in a two-dimensional turbulent channel flow[END_REF], Tiselj et al. [START_REF] Tiselj | Effect of wall boundary condition on scalar transfer in a fully developed turbulent flume[END_REF]).

The present DNS is performed with a time step of 0.008 ν/u 2 τ and one snapshot is taken every 1000 time steps. 700 snapshots are used to reconstruct the statistics, corresponding to an averaging time of 5600 ν/u 2 τ .

DNS is widely used to produce reference data. However, as pointed out by Oliver et al. [START_REF] Oliver | Estimating uncertainties in statistics computed from direct numerical simulation[END_REF], statistics obtained from DNS contain non-trivial errors.

Errors arise mainly from the discretization of the equations and from the finite statistical sampling. As our code is based on a pseudo-spectral method and we use a fine grid, the spatial discretization error of the Navier-Stokes equation is rather low. Thus, the main source of error in our statistics is the sampling error, which we estimate with the open-source library ar (https://github. com/RhysU/ar, see Oliver et al. [START_REF] Oliver | Estimating uncertainties in statistics computed from direct numerical simulation[END_REF]).

The program arsel, which ships with the library, is used as follows. For each statistical quantity of interest, the instantaneous value is extracted from each snapshot and averaged in the homogeneous directions. The resulting time series are used to fit autoregressive models with arsel using the --subtract-mean option. The output of arsel is then processed to extract the average profile and the sampling error profile. Herein, the sampling error is given as 2 standard deviations. Given the domain extension, the number of modes used and the time step, the probability that the real average and its estimation differ for less than 2 standard deviations is 95%.

As this method can be applied for any statistical quantity, we are able to estimate the sampling error for the mean temperature, the temperature dissipation rate or its variance. For a given variable φ, the associated sampling error is e φ . We also define the dimensionless normalized sampling error E [φ] by

E [φ] (y) = e φ (y) e φ (y) + abs φ (y) (5)
This normalized sampling error is bounded in [0, 1], the upper bound being reached when the averaged value abs φ (y) is exactly zero. As most of the time the averaged value is large compared to the sampling error, E [φ] measures the relative uncertainty.

Preliminary results on coarser grids

The sampling error on u 2 τ , whose theoretical value is 1 herein, and ∂ y θ (y = ±1) is illustrated with the four DNS studies on coarser grids described in the Table 2. The results in the large computational box combined with the long and short averaging time are taken from the analyses of Tiselj and Cizelj [START_REF] Tiselj | Dns of turbulent channel flow with conjugate heat transfer at prandtl number 0.01[END_REF], which were performed at Re τ = 180, P r = 1 and with slightly different thermal boundary conditions: both channel walls were heated and a forcing term was present in the passive scalar equation (4). The small box results are part of the present analysis: a separate run over longer time interval was made only with a resolution, which corresponds to the case 1KMM in Table 1.

Label Re τ L + x L + y L + z ∆ + x ∆ + y ,max ∆ + z Averaging time in wall units u 2 τ ± e u 2 τ ∂ y θ ± e ∂yθ 18L 18S
The first impression about the statistical uncertainty can be obtained from the two rightmost columns of Table 2 where the squared friction velocities and temperature derivatives at both channel walls are shown. As one might expect, deviation from the theoretical friction velocity of 1 is larger in the smaller computational domain and at shorter averaging time intervals. As for the temperature derivative, shorter averaging time also leads to higher sampling error.

Regarding cases Re τ = 160, we observe a relatively large difference on the temperature derivative at both walls. Although this difference gets lower when increasing the duration of the simulation, it is probably specific to smaller computational domains, as reported by Jimenez and Moin [START_REF] Jiménez | The minimal flow unit in near-wall turbulence[END_REF] for velocity statistics.

The case 16S in Table 2 is comparable to the case 1KMM in Table 1. Both cases correspond to an averaging time roughly 2 times longer than the one in Galantucci and Quadrio [START_REF] Galantucci | Very fine near-wall structures in turbulent scalar mixing[END_REF], with snapshots taken at a similar frequency for all those cases. Thus, the case 16S gives a lower bound for the sampling error in Galantucci and Quadrio [START_REF] Galantucci | Very fine near-wall structures in turbulent scalar mixing[END_REF]. The Figure 1 shows the normalized sampling error associated with various statistical quantities for the cases present in Table 2. The clear trend in Figure 1 is that higher order moments are associated with higher sampling errors. Regarding the plots on the middle row, the normalized sampling error for the temperature variance and for the mean dissipation rate have a similar amplitude and involve turbulent quantities to the power 2. Qualitatively, the situation is similar on the plots in the bottom row, which involve turbulent quantities to the power 4.

Regarding case 16S, the magnitude of the normalized sampling error is around 5% for E [θ], except in the middle of the channel, where it gets close to 100% as the average temperature is zero. It also drops down to zero at the walls as the temperature is imposed there (Dirichlet boundary condition, θ = ±1).

For cases 18S and 18L, the imposed temperature at the walls is zero, so both the sampling error and the average temperature tend to zero there. The resulting normalized sampling error E [θ] at the wall for those cases is around 0.2%.

Focusing on case 16S for higher order integral statistics, the normalized sampling error reaches almost 10% for E θ 2 and E [ε θ ], and even 15% for E θ 4 and E ε 2 θ . These results show that the differences of 7% on the average of ε θ on different grid reported in Galantucci and Quadrio [START_REF] Galantucci | Very fine near-wall structures in turbulent scalar mixing[END_REF] is well within the sampling error as their simulations have the same domain extension but a duration halved compared with the case 16S.

Using the Central Limit Theorem applied to weakly correlated samples described in Oliver et al. [START_REF] Oliver | Estimating uncertainties in statistics computed from direct numerical simulation[END_REF] (Theorem 27.4 in Billingsley [START_REF] Billingsley | Probability and measure[END_REF]), the sampling error decreases with the square root of the number of samples. For a given sampling frequency, the error is therefore halved when the duration of the simulation (and the number of samples) is increased by a factor 4. This holds quite well when comparing cases 16S and 16L, or 18S and 18L. Overall, significant reduction of the sampling error can be achieved using extended computational domains and / or simulations with a longer duration.

θ θ 2 θ 3 θ 4 ε θ ε 2 θ 1KMM 40 *

Impact of unresolved smaller scales on statistics

In this section, we investigate the impact of unresolved smaller scales on various quantities. First, we investigate dimensionless bulk quantities. Using the scalar with the finest resolution (6S) as a reference, we can compute a quadratic relative error. For instance, the quadratic relative error for the average temperature on the grid 1KMM is defined by:

y=1 y=-1 θ 1KM M (r) -θ 6S (r) 2 dr y=1 y=-1 θ 6S (r) 2 dr (6)
In Table 3, the quadratic relative error is given for various quantities. Looking at the relative error for the moments of the temperature (θ, θ 2 , θ 3 and θ 4 ), one can extract two trends. First, for all resolutions considered, the higher the moment, the higher the relative error. Second, the finer the resolution, the closer to the case 6S. The resolution 1KMM seems sufficient to accurately estimate the average temperature. However, accurate estimation of higher order moments requires a finer resolution as unresolved smaller scales have a higher impact on them. The situation is similar for ε θ and ε 2 θ : the first is relatively well estimated by the coarsest grid while the second has a stronger dependence on unresolved smaller scales.

It is important to stress that the discretization error in Table 3 is typically several decades lower than the sampling error. However, all the scalar fields are transported by the same velocity field. Therefore, they are strongly correlated and are all plagued by the same sampling error. This allows us to exhibit the impact of unresolved smaller scales on statistics, even though this impact is relatively small compared with the sampling error.

For the turbulent channel flow studied here, at the wall, the temperature is imposed (θ = ±1) and the temperature evolution equation 4 contains no source term. As a result, the exact average temperature derivative at the wall is plagued by statistical uncertainty and depends on the spatial resolution. We define the friction temperature θ τ and the reference temperature θ ref with:

θ τ = 1 Re τ P r ∂ y θ ± L y 2 and θ ref = 1 Re τ P r θ (L y /2) -θ (-L y /2) L y (7)
The theoretical value of the friction velocity, which is equal to one here, produces a θ ref not plagued by statistical uncertainty, while θ τ is. Therefore, we use θ ref to rescale the temperature statistics. Spatial derivatives are expressed in wall-units using the theoretical value of the friction velocity. For convenience, plots are using a wall-normal distance in computational units (y ∈ [-1, 1]).

Complementary with the bulk quantities in Table 3 and the integral profiles of ε θ or its variance, we investigate the anisotropy of the fluctuating temperature gradient. This is done using the barycentric map defined by Banerjee et al. [START_REF] Banerjee | Presentation of anisotropy properties of turbulence, invariants versus eigenvalue approaches[END_REF].

Given a symmetric positive semi-definite tensor M ij , the associated (traceless)

anisotropy tensor is

M ij = M ij M kk - 1 3 δ i,j (8)
The eigenvalues of M ij , λ k , verify λ 1 ≥ λ 2 ≥ λ 3 . Based on those eigenvalues, Banerjee et al. [START_REF] Banerjee | Presentation of anisotropy properties of turbulence, invariants versus eigenvalue approaches[END_REF] define 3 limiting states forming a basis: one-component (1C), two-component (2C) and three-component (3C). The number of components corresponds to the number of non-zero eigenvalues of M ij . In the state 3C, all three eigenvalues are equal: the tensor is isotropic. In the state 2C, one eigenvalue is negligible while the two others are equal: the anisotropy is twodimensional. In the state 1C, two eigenvalues are negligible: the anisotropy is one-dimensional.

In the barycentric map, each limiting state is the vertex of an equilateral triangle. For any symmetric positive semi-definite tensor M ij , there is a unique point inside this equilateral triangle with coordinates

[x, y] = (λ 1 -λ 2 ) [x 1C , y 1C ] + 2 (λ 2 -λ 3 ) [x 2C , y 2C ] + (3λ 3 + 1) [x 3C , y 3C ] (9)
where [x iC , y iC ] is the coordinate of the vertex corresponding to the limiting state iC. Following Lumley and Newman [START_REF] Lumley | The return to isotropy of homogeneous turbulence[END_REF], it is possible to visualize the anisotropy of a traceless tensor using its second and third invariants. However, the resulting invariant map is distorted as those invariants are nonlinear functions of the tensors elements. The barycentric map is connected to the invariant map with a bijection and does not suffer from this bias. This is the main motivation for the present choice; interested reader will find a deeper analysis of the barycentric map in Banerjee et al. [START_REF] Banerjee | Presentation of anisotropy properties of turbulence, invariants versus eigenvalue approaches[END_REF].

The symmetric positive semi-definite tensors investigated here are related to the temperature derivatives. The first tensor investigated, ∂θ

= ∂ i θ ∂ j θ ,
is relevant to this study as its trace is proportional to the scalar dissipation rate ε θ and its anisotropy is related to the anisotropy of the fluctuating temperature gradient. The other tensor investigated is the variance of ∂θ :

σ ∂θ = [∂ i θ ∂ j θ ] 2 -∂ i θ ∂ j θ 2 .
Locations where the fluctuating temperature gradient is isotropic correspond to points of the barycentric map closer to the vertex 3C of the triangle. The further away from this vertex, the higher the anisotropy.

Points closer to 1C or 2C correspond to one-dimensional or two-dimensional fluctuations, respectively. It is important to stress that the barycentric map does not provide a complete description of the anisotropy as it obliterates the eigenvectors associated with the eigenvalues.

The left frame of Figure 2 is in agreement with Table 3; the profiles of ε θ do not show any significant impact of unresolved smaller scales in 1KMM and 3S results. The situation is similar for the barycentric map of the anisotropy of ∂θ (right frame of Figure 2). Due to the imposed temperature at the wall, the anisotropy at the boundary has only one component (limiting state 1C). Then, the closer to the middle of the channel, the closer to isotropy (limiting state 3C).

Focusing on the left frame of Figure 2, one may notice that the profiles of ε θ are not exactly symetric. This is not an issue here as this is well within the sampling error, which is around 10% for ε θ , according to Figure 1. However, performing a simulation with a duration long enough to have a negligible sam-260 pling error, one may obtain a profile that still lacks symetry. This is specific to the combination of smaller domains and periodicity used here: global mixing between the lower and upper half of the domain is not perfect due to the absence of large scale structures. The left frame of Figure 3 is also in agreement with Table 3; the variance of ε θ is significantly impacted by unresolved smaller scales, and the coarsest grid underestimates this quantity, especially outside of the viscous sublayer where the discretization error is as large as the sampling error. Closer to the wall, this impact is less visible and all spatial resolutions produce a similar variance of ε θ . The barycentric map of the anisotropy of σ ∂θ also shows the impact of unresolved smaller scales. In the middle of the channel (y = 0), the coarsest grid gives a similar weighting for all limiting states (1C, 2C and 3C), this correspond to the part of 1KMM closer to the middle of the triangle. At this location of the channel (y = 0), finer grids give a similar weighting only for the limiting states 1C and 3C, but a negligible weight to state 2C, this correspond to the part of 6S and 3S closer to the edge connecting 1C and 3C. Integral quantities provide only a limited amount of information and spectra in Figure 4 are more insightful. They are taken at the middle of the channel (y = 0). The spanwise spectra show that the spanwise contribution to ε θ is well resolved even on the coarsest grid. However, this grid is not isotropic and it has extended cells in the streamwise direction (∆ + x = 17, ∆ + z = 6 for 1KMM). This leads to a pronounced bump on the streamwise spectrum of 1KMM at the highest resolved wavenumbers (k + x ≈ 0.2). As one can notice, the streamwise spectrum of 3S also has a bump at the highest resolved frequency (k + x ≈ 0.55). However, regarding the latter case, this has a negligible impact on the statistics of ε θ because the bumps correspond to very high frequencies for which the streamwise contribution to the dissipation is two decades smaller than its peak value around k +

x ≈ 0.05. This does not hold for case 1KMM: the bump on the streamwise spectrum is present at lower frequencies, for which the streamwise contribution to the dissipation has a magnitude similar to the one at the peak. This explains the severe impact of unresolved smaller scales on the variance of ε θ exhibited in Figure 3.

Conclusion

The present investigation demonstrates the impact of a limited resolution of the smaller scales on the scalar dissipation rate, its variance and on the anisotropy of the fluctuating scalar gradient. In a preliminary part, we showed that simulations of the turbulent channel flow based on smaller domains and with a limited duration (often used to perform grid sensitivity analysis) are plagued by a relatively large sampling error, especially for higher order moments like the scalar dissipation rate or its variance. This sampling error is often larger than the discretization error, thus preventing accurate estimation of the latter, except for simulations with an extended domain and duration.

In the present study, this situation is overcomed by using a single velocity field to transport several passive scalars with different resolutions: all the scalars are strongly coupled by the velocity field and plagued by the same sampling error. The resulting difference between the scalars is solely due to the discretization error, which can thus be estimated. Another workaround could be to perform simulations with a duration long enough so that the sampling error becomes negligible compared with the discretization error, or to perform many simulations with slightly distinct initial conditions, followed by ensemble averaging.

Regarding the key integral turbulence statistics of the passive scalar investigated here, our results confirm that previous DNS based on a pseudo-spectral method with the resolution 1KMM (Kim et. al. [2], Kasagi et al. [START_REF] Kasagi | Direct numerical simulation of the passive scalar field in a two-dimensional turbulent channel flow[END_REF], Tiselj et al. [START_REF] Tiselj | Effect of wall boundary condition on scalar transfer in a fully developed turbulent flume[END_REF]) can estimate the average temperature, its variance and the average scalar dissipation rate with a relatively low discretization error, as indicated Table 3. However, they can only provide qualitative estimation for higher order moments. To accurately estimate those higher order moments, finer grids must be used. Regarding the variance of the scalar dissipation rate, or the anisotropy of the fluctuating scalar gradient, the grid recommended by Vreman and Kuerten [START_REF] Vreman | Comparison of direct numerical simulation databases of turbulent channel flow at re τ = 180[END_REF] is sufficient for pseudo-spectral schemes (6 th order compact finite-difference schemes typically require a grid 4/3 finer). Our results do not preclude existence of scalar structures or events, which are too small to be resolved with the Vreman and Kuerten resolution. However, our results show that these events might significantly affect only the integral turbulence statistics higher than fourth-order. It is important to stress that a small event can impact part of the PDF of a turbulent quantity without significantly affecting the associated integral value.

Our results also suggest that unresolved smaller scales have a stronger impact outside of the viscous sublayer, and have almost no impact closer to the wall. Strictly speaking, this observation is valid only for the Dirichlet boundary condition investigated here (imposed temperature), and may not hold for a Neumann boundary condition (imposed heat flux), a Robin boundary condition (heat exchange coefficient) or conjugate heat-transfer (fluid-solid thermal coupling). Regarding the latter case, the present analysis will provide a solid ground for future investigations on the discontinuity of ε θ at the fluid-solid interface.

- 1 .Figure 1 :

 11 Figure 1: Normalized sampling error for the cases in Table 2. Top: E [θ]. Middle: E θ 2 (left) and E [ε θ ] (right). Bottom: E θ 4 (left) and E ε 2 θ (right).

Figure 2 :

 2 Figure 2: Dissipation rate ε θ (left) and barycentric map associated with ∂θ (right).

Figure 3 :

 3 Figure 3: 95% confidence interval of the variance of the dissipation rate ε θ (left) and barycentric map associated with σ ∂θ (right).

Figure 4 :

 4 Figure 4: Spectra of the streamwise (left) and spanwise (right) contributions to the scalar dissipation rate ε θ at y = 0 (middle of the channel).

Figure 5 :

 5 Figure 5: Instantaneous scalar gradient ( ∂ i θ 2 ) in computational units at y = 0 (middle of the channel). Logarithmic colormap. Left: 1KMM. Right: 6S.

Figure 5

 5 Figure 5 illustrates the impact of the under-resolution on the gradient of the scalar. On the one hand, both fields exhibit similar structures at larger scales, in agreement with the ability of the grid 1KMM to estimate the average scalar dissipation rate, as shown in Figure 2. On the other hand, smaller scales are visibly misrepresented on the coarser grid, as the picture on the left is plagued by the oscillations, which correspond to the wavenumbers of around k + = 0.2 that are overpredicted in 1KMM results shown in the left frame of

Figure 4 .

 4 Figure 4. This is consistent with the previous results: the grid 1KMM is too coarse to properly resolve the sharp gradients which separate dissipative and non-dissipative regions.

Table 1 :

 1 Spatial resolution for the transported passive scalar fields.

	19,

Table 2 :

 2 Sampling error on friction velocity and temperature derivative at the boundaries.

Table 3 :

 3 log 10 of the quadratic relative error for various quantities.

		129 * 60	-7.18 -5.02 -4.42 -3.67 -6.68 -1.26
	2GQC	64 * 129 * 64	-9.15 -7.14 -6.50 -5.43 -7.82 -1.91
	3S	112 * 129 * 80 -11.8 -10.5 -8.57 -8.11 -9.54 -3.25
	4VK	112 * 181 * 80 -11.8 -10.5 -8.58 -8.13 -9.98 -3.29
	5GQM	360 * 129 * 180 -13.9 -12.0 -10.0 -10.8 -9.51 -5.53

Acknowledgements

This work was financially supported by the research project of the Slovenian Research Agency P2-0026 and by the EDF-JSI collaboration, project PR-07184.

Data associated with the present paper are available online at http://dx. doi.org/10.17632/mn74gv69wn.1 and https://repo.ijs.si/CFLAG/sml-scl under the GNU GPL v3 licence.