Bruno Chaouat

Extension of the Partially Integrated Transport Modeling method to the simulation of passive scalar turbulent fluctuations at various Prandtl numbers

The present work aims to extend the hybrid non zonal RANS/LES partially integrated transport modeling (PITM) method to turbulent flows in the presence of passive scalar contaminant for simulating large scales of turbulent flows. Focussing on the methodological aspects, we derive the basic transport equations both for the scalar variance of fluctuations and dissipation-rate of the variance. The basis of the method was introduced in references [R. Schiestel and A. Dejoan, "Towards a new partially integrated transport model for coarse grid and unsteady turbulent flow simulations ", Theor. Comput. Fluid Dyn. 18, 443 (2005)] and [B. Chaouat and R. Schiestel, "A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows ", Phys. Fluids 17, 065106 (2005)]. It provides a continuous approach for hybrid Reynolds averaged Navier-Stokes equations-large eddy simulation (RANS-LES) with seamless coupling between RANS and LES regions. The main motivation is to simulate accurately in LES mode performed on coarse grids, scalar fluctuation fields, in addition to mean scalar fields.

As known, the knowledge of the rms scalar fluctuations is often involved in handling practical engineering and geophysical flows. Like in dynamical equations, it is found that the coefficient appearing in the destruction term of the dissipation-rate of the scalar variance is a function of the cutoff-wave number and also of the Prandtl number and the Reynolds number. Depending on the value of the Prandtl-number, different expressions of this function have been derived according to the relevant physics in the wave number space. Finally, numerical simulations of fully turbulent flows including passive scalar transport fields have been performed on several meshes of medium and coarse grid resolutions at the Reynolds number R τ = 395 for the Prandtl numbers P r = 0.1, 1 and 10, respectively associated with heat transfer of liquid metals, gas and water for illustrating the theoretical development made on PITM. As expected, the PITM method provides satisfactory results in good agreement with data of direct numerical simulations. From a general point of view, this work opens new routes of modeling and simulation of turbulent flows including a passive scalar with a drastic reduction of the computational time and memory in term of number of grid points in comparison with the demanding resources of highly resolved LES.

Introduction

Mathematical turbulence modeling methods such as Reynolds averaged Navier-Stokes (RANS) or large eddy simulations (LES) methods have been proposed independently from each other for simulating turbulent flows. Each of them has their respective advantage and limitation and their specific field of application so that they should be considered as complementary tools in computational fluid dynamics.

In that sense, the most appropriate tool does not depend on the intrinsic performances of the method itself but more precisely on both the required computational ressource and nature of the flow. RANS models have been developed initially for simulating flows that evolve slowly in time with a reasonable computational cost and are rather used in industry for determining practical data such as the pressure and heat coefficients used in the optimisation of engineering processes while LES demanding large computational cost is rather used for simulating unsteady flows subjected to high frequencies and are often studied in research laboratories with emphasis on fundamental aspects and the turbulence structures. Since the past two-decade, hybrid RANS-LES have been proposed to take benefit of both RANS and LES methods [START_REF] Fröhlich | Hybrid LES/RANS methods for the simulation of turbulent flows[END_REF]Von Terzi, 2008, Chaouat, 2017). Among these numerous hybrid RANS/LES methods, the partially integrated transport modeling (PITM) is an advanced method that has been developed in previous publications [START_REF] Schiestel | Towards a new partially integrated transport model for coarse grid and unsteady turbulent flow simulations[END_REF][START_REF] Chaouat | A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows[END_REF][START_REF] Gatski | Second-moment and scalar flux representations in engineering and geophysical flows[END_REF][START_REF] Chaouat | Simulation of turbulent rotating flows using a subfilter scale stress model derived from the partially integrated transport modeling method[END_REF], 2013;[START_REF] Chaouat | Subfilter scale transport model for hybrid RANS/LES simulations applied to a complex bounded flow[END_REF][START_REF] Chaouat | Simulation of turbulent rotating flows using a subfilter scale stress model derived from the partially integrated transport modeling method[END_REF]Chaouat, , 2017b) ) for simulating large eddy scale of turbulent flows out of spectral equilibrium performed on coarse grids with a mesh step corresponding to a spectral cutoff wave number that may be located before the inertial zone of the Kolmogorov range. This type of approach allows continuous hybrid non-zonal simulations between near RANS regions and LES regions with seamless coupling. The PITM model reduces to RANS when the grid step increases, if the flow is homogeneous in the three directions. Depending on the value of the spectral splitting parameter, usually linked to the numerical mesh, each part of the fluctuating turbulent flow is partially modeled and partially simulated like in the LES method. There is a unique link between the filter width or the splitting wavenumber and the turbulent energy ratio of subfilter energy to simulated energy. The PITM models were the first to introduce this practice [START_REF] Schiestel | Towards a new partially integrated transport model for coarse grid and unsteady turbulent flow simulations[END_REF][START_REF] Chaouat | A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows[END_REF] and some other models of the literature, like the partially averaged Navier-Stokes (PANS) models recently borrowed this formulation in order to improve their formalism [START_REF] Foroutan | A partially averaged Navier Stokes model for the simulation of turbulent swirling flow with vortex breakdown[END_REF]. The method is originally developed on a spectral basis, but recently a variational analysis has been conducted to handle inhomogeneous flows [START_REF] Heinz | The large eddy simulation capability of Reynolds-averaged Navier-Stokes equations: Analytical results[END_REF]. The main feature of this method lies in its general character that allows to convert almost any usual one point RANS transport closure into a corresponding hybrid non-zonal subfilter-scale model for performing large eddy scale simulations. So, PITM can be applied to two-equation turbulence models but also to more advanced turbulence models such as Reynolds stress models developed in the framework of second moment closure (SMC) [START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF][START_REF] Hanjalic | Modelling turbulence in engineering and the environment. Second-moment route to closure[END_REF].

The present work deals with the extension of the PITM method to passive scalar fluctuations. We shall give emphasis in the theoretical development of the methodology and derivation of equations applied to scalar variance from a physical standpoint. Applications will be restricted here to the case of the plane channel flow with associated heat transfer in order to demonstrate the satisfactory behavior of the model. More extensive and varied applications will be considered in subsequent further works.

This work is highly motivated by the fact that passive scalar fields are often encountered both in nature including geophysical and environmental flows as well as engineering flows in industrial processes. In RANS and LES simulations involving scalar fields for temperature or concentration, it is of primary importance to accurately model the fluxes of the scalar-velocity correlations that play a major role in the determination in heat and mass transfer in turbulent flows. Indeed the mean scalar field equations are directly coupled with the turbulence dynamic equations. But, the knowledge of the rms fluctuations of the scalar field may also be of importance in several practical cases when the high scalar fluctuations can impact the wall structures. Structures exposed to such temperature fluctuations may suffer from thermal fatigue damage. This problem may appear in the mixing of two flows at different temperatures merging at a T-junction. For instance, among others, thermal striping can arise in certain liquid metal-cooled reactors. Because this may result in thermal fatigue cracking, the problem deserves study.

In the present work, we shall give a special focus to the study of the variance of fluctuations of the scalar, keeping in mind that scalar transfer problems are inherently interrelated. The same general form for the turbulent scalar equations invites us to consider any kind of scalar. Overall, among the large variety of passive scalar, we can for instance point out the temperature when it is possible to neglect the dynamic effect of density variations but in a more general way, it can be also any passive contaminant concentration transported by the flow itself that deserves physical investigations.

After briefly recalling the basics of the PITM method, we will extend it to the study of scalar contaminant fluctuations. In the present case according to previous works of the authors [START_REF] Chaouat | A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows[END_REF][START_REF] Gatski | Second-moment and scalar flux representations in engineering and geophysical flows[END_REF][START_REF] Chaouat | Simulation of turbulent rotating flows using a subfilter scale stress model derived from the partially integrated transport modeling method[END_REF], 2013) and [START_REF] Chaouat | Subfilter scale transport model for hybrid RANS/LES simulations applied to a complex bounded flow[END_REF][START_REF] Chaouat | Simulation of turbulent rotating flows using a subfilter scale stress model derived from the partially integrated transport modeling method[END_REF]Chaouat ( , 2017)), we will retain the second moment closure for modeling the subfilter scale stresses considering that this is the appropriate level of closure to get accurate results of complex turbulent flow fields. However a first order closure based on diffusivity models will be sufficient in a first step for modeling the scalar fluxes involved in the production of variance equation for the subfilter passive scalar [START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF][START_REF] Hanjalic | Modelling turbulence in engineering and the environment. Second-moment route to closure[END_REF]. There are several arguments in favor of this choice. Firstly, the turbulent velocity field is here decoupled from the passive scalar field. Secondly, tractable equations can be solved easily from an engineering point of view in both research and engineering CFD codes. Hence, this strategy to choose different levels of closure for each type of flow field is justified here both from physical and numerical points of view while allowing fair results for the scalar field. In this work, we will give emphasis on the PITM extension to scalar transport that comes from the new theoretical developments in the dissipation rate of variance equation. As this equation is common both to first order and higher order closure, it would be a simple matter to extend the present work to second order passive scalar closure.

2

The basics of the PITM method 2.1 General framework : the levels of description

The PITM method is funded originally on the spectral equations for quasi-homogeneous turbulence.

With an hypothesis of tangent homogeneous space [START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF] it becomes then possible to deal with non-homogeneous flows with a good approximation for practice. The steps used in the development of the model introduce several level of description which are made clear in the synoptic (Table 1) for dynamic and scalar turbulent quantities. In shear driven turbulent flows in the presence of a passive scalar, the transport equation for the turbulent stresses and the turbulent scalar fluxes Φ ij and Φ iθ , respectively, appearing in lines 3 and 4 of table 1 are the most important tensors in second order closures to describe turbulent transfer. These equations, at all level of description, contain important redistribution terms linked to pressure fluctuations that have been a challenge in spectral closures and in RANS closures as well. In the turbulent stresses closure, the redistribution term of turbulence energy is usually modeled using several contributions including a linear contribution (rapid part involving the mean velocity gradients) and a non-linear contribution (slow part or Rotta term) along with a wall contribution. In the equation for the scalar fluxes, a similar approach can be used, introducing the corresponding three contributions. The study of these term in the spectral context has been extensively developed [START_REF] Cambon | Spectral modelling of homogeneous non-isotropic turbulence[END_REF]. If the linear terms are closed at level 1, a parametrization of angular dependence with the wave vector κ is necessary at level 2. As a consequence of spherical averaging, a closure for linear terms is required like in RANS models. A more advanced representation of the linear terms is worked out [START_REF] Mons | A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically averaged descriptors[END_REF]) by introducing two anisotropy tensors for the polarization anisotropy and for the directional anisotropy in wave vector, while the non-linear transfer is obtained from eddy damped quasi-normal Markovian (EDQNM) hypothesis. This spectral modeling approach has been extended also to passive scalar transport [START_REF] Briard | Spectral modelling for passive scalar dynamics in homogeneous anisotropic turbulence[END_REF].

In this model, the equations for temperature variance as well as the turbulent thermal fluxes have been closed using directional and polarization parameters in homogeneous anisotropic fields. For the subfilter equations (level 3) the pressure-strain closures are inspired from RANS models and will be used in this paper both for the slow and rapid parts. So, in the present simulations, the subfilter stress tensor (τ ij ) sf s is numerically calculated from its transport equation considering the turbulence

E(X, κ, t) E(X, κ, t) k sf s (X, t) k sf s (X)
Variance of passive scalar fluctuations

E θ (X, κ, t) E θ (X, κ, t) k θsf s (X, t) k θsf s (X) Double velocity corre- lations Φ ij (X, κ, t) ϕ ij (X, κ, t) (τ ij ) sf s (X, t) (τ ij ) sf s (X) Velocity passive scalar correlations Φ iθ (X, κ, t) ϕ iθ (X, κ, t) (τ iθ ) sf s (X, t) (τ iθ ) sf s (X)
spectral tensor ϕ ij . The same type of terms appear in the equation for the subfilter scalar fluxes (τ iθ ) sf s that will be considered subsequently using the turbulence scalar spectral vector fluxes ϕ iθ . It is worth noting that in passive scalar transport, the scalar variance k θ as well as its dissipation rate ǫ θ are not at all involved in the dynamic equations. But the reverse is not true, the equations for these quantities are dependent on many other turbulent correlations including the scalar-velocity correlations in addition to the dynamical correlations. In particular, the production term in the scalar variance equation involves the aforementioned scalar-velocity correlation tensor corresponding physically to the turbulent scalar fluxes ϕ iθ . A modeled transport equation for these fluxes is of utmost importance in the study of turbulent scalar transport and would be required here to stick to the second order of closure.

However, considering that our present focus is mainly devoted to the study of scalar fluctuations, we shall postpone this requirement to a further study. We will retain here a gradient approximation of this term using a tensorial diffusivity. The main concern of the present paper is therefore to extend the spectral methods already developed by the present authors in dynamic turbulence to the fluctuations of the transported scalar.

Spectrum splitting and partial integration

In a general way, each flow variable φ can be decomposed into a statistical mean value φ and a fluctuating turbulent part φ ′ which is subdivided itself into two ranks of fluctuating parts, a large scale part and a small scale part (Schiestel, 1983a(Schiestel, , 1983b[START_REF] Schiestel | Multiple-time scale modeling of turbulent flows in one point closures[END_REF][START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF]. This decomposition can be written as φ = φ + φ < + φ > where φ < and φ > can be defined by the spectral operators

φ < (ξ) = 0<|κ|<κc φ ′ (κ) exp (jκξ)dκ, (1)
and

φ > (ξ) = κc<|κ|<∞ φ ′ (κ) exp (jκξ)dκ, (2)
in which κ c , here, denotes the first splitting wave number or spectral cutoff. From a physical standpoint, the PITM method finds its basic foundation in the spectral space of wave numbers considering the production, transfer and dissipation processes of energy acting in each spectral wave number range of the spectrum. Different steps must be undertaken to get the basic equations that constitute the cornerstone of the PITM method. In the first step, for any turbulent variable, we develop the transport equation for the two-point fluctuating correlations in the physical space using the concept of the tangent homogeneous space (Chaouat andSchiestel, 2007, Chaouat, 2017) accounting for the midway position X = (x A + x B )/2 and the vector difference ξ = x Bx A between the two points A and B. In the second step, we take the Fourier transform of this equation. In the third step, we apply the spherical mean operator [START_REF] Cambon | Spectral modelling of homogeneous non-isotropic turbulence[END_REF][START_REF] Schiestel | Multiple-time scale modeling of turbulent flows in one point closures[END_REF][START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF] to get an equation that depends only on the scalar wave number and not anymore on the vector wave number itself. In the last step, we finally return to the physical space from partial integration of the spectrum in the wave number ranges [κ c , κ d ] to get the corresponding transport equation for the correlation of the fluctuating variable φ > φ > in the physical space. This mathematical physics formalism will be followed in the present work to derive equations for the subfilter scale of the passive scalar variance, as already performed for the dynamics of the flow field [START_REF] Schiestel | Multiple-time scale modeling of turbulent flows in one point closures[END_REF][START_REF] Schiestel | Towards a new partially integrated transport model for coarse grid and unsteady turbulent flow simulations[END_REF][START_REF] Chaouat | A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows[END_REF][START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF][START_REF] Gatski | Second-moment and scalar flux representations in engineering and geophysical flows[END_REF]) by using the shape of spectra as studied by [START_REF] Batchelor | Small scale variation of convected quantities like temperature in a turbulent fluid. Part 1. General discussion and the case of small conductivity[END_REF] and reinterpreted recently by [START_REF] Warhaft | Passive scalars in turbulent flows[END_REF].

2.3

Recalling of turbulent velocity field equations

This section recalls the main concepts of the PITM method that allows to derive the major governing equations for the velocity field. Indeed this is not useless because the extension to scalar transport has to rely upon the dynamic equations. The general approach is founded on the technique of partial spectral integration. First consider splitting the turbulence energy spectrum as sketched in Fig. 1 by introducing a cut-off wavenumber. It can be viewed as a filtering operator in spectral space. We suppose for the present purpose that turbulence is fully developed at high Reynolds number and we will consider low Reynolds number modifications later on. As usually, the Reynolds stress tensor is defined

LogE k k k c d k k F c
τ ij (X, t) = u ′ i u ′ j (X, t) (3)
where u ′ i denotes the usual fluctuating velocity appearing in the statistical treatment and X denotes the midway-position between the two points and t is the time. The starting point is the transport equation of the spherical mean of the Fourier transform of the two-point correlation tensor of the fluctuating velocities denoted ϕ ij (X, κ, t) = u ′ i u ′ j (X) ∆ (κ, t) as follows [START_REF] Hinze | Turbulence[END_REF][START_REF] Schiestel | Multiple-time scale modeling of turbulent flows in one point closures[END_REF][START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF])

∂ϕ ij (X, κ, t) ∂t + u j (X) ∂ϕ ij (X, κ, t) ∂X j = P ij (X, κ, t) +T ij (X, κ, t) + Ψ ij (X, κ, t) + J ij (X, κ, t) -E ij (X, κ, t) (4)
where P ij , T ij , Ψ ij , J ij , and E ij are respectively, the production, transfer, redistribution, diffusion and dissipation terms. All the expressions of the terms appearing in Eq. ( 4) can be identified using a mathematical-physics formalism developed in the spectral space [START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF]. The transport equation for the turbulence energy spectrum E = ϕ jj /2 is simply obtained by tensorial contraction of Eq. ( 4) leading to

∂E(X, κ, t) ∂t + u j (X) ∂E ij (X, κ, t) ∂X j = P(X, κ, t) +T (X, κ, t) + J (X, κ, t) -E(X, κ, t) (5)
where P = P mm /2, T = T mm /2, J = J mm /2 and finally, E = E mm /2. Partial integration can be performed on wave number ranges such as [0, κ c ], [κ c , κ d ] and [κ d , ∞[ where κ c denotes the cutoff wave number given by the grid size ∆ and κ d is the dissipative wave number located at the far end of the inertial range of the spectrum assuming that the energy pertaining to higher wave numbers is negligible.

Taking into account the significant processes of turbulence acting in each spectral zone of the density spectrum, it is then possible with some calculus to derive the equations associated with the subfilterscale and dissipation rate as demonstrated step by step by [START_REF] Schiestel | Towards a new partially integrated transport model for coarse grid and unsteady turbulent flow simulations[END_REF] as well as Chaouat andSchiestel (2005, 2012).

Equation set for the turbulent dynamic field

∂(τ ij ) sf s ∂t + ∂ ∂x k (ū k (τ ij ) sf s ) = (P ij ) sf s + (Π ij ) sf s + J ij -ǫ ij (6)
where the terms appearing in the right-hand side of this equation are identified as subfilter production, redistribution, diffusion and dissipation, respectively. The transport equation for the subfilter turbulent energy is obtained as half the trace of Eq. ( 6)

∂k sf s ∂t + ∂ ∂x k (ū k k sf s ) = P sf s + J -ǫ (7)
where P sf s = (P mm ) sf s /2, J = J mm /2, ǫ = ǫ mm /2. The production term (P ij ) sf s accounts for the interaction between the subfilter stresses and the filtered velocity gradients

(P 1 ij ) sf s = -(τ ik ) sf s ∂ ūj ∂x k -(τ jk ) sf s ∂ ūi ∂x k (8)
The redistribution term (Π ij ) sf s appearing in Eq. ( 6) is decomposed into a slow part (Π 1 ij ) sf s that characterizes the return to isotropy due to the action of subgrid turbulence on itself

(Π 1 ij ) sf s = -c 1sf s ǫ k sf s (τ ij ) sf s - 1 3 (τ mm ) sf s δ ij (9)
and a rapid part, (Π 2 ij ) sf s that describes the action of the filtered velocity gradients

(Π 2 ij ) sf s = -c 2 (P ij ) sf s - 1 3 (P mm ) sf s δ ij (10)
where c 1sf s = c 1 α(η) is an increasing function of the parameter η to strengthen the return to isotropy for large wave numbers, c 1 and c 2 are constant coefficients. The closure of the rapid term in the pressurestrain correlations is inspired from the stress transport model developed by [START_REF] Launder | Progress in the development of a Reynolds stress turbulence closure[END_REF] referring to the so-called "isotropization of production "which is a simplified form of the "quasi isotropic "version proposed in the same paper. A similar modeling was obtained by Jeandel et al. [START_REF] Jeandel | Modeling methods in physical and spectral spaces[END_REF] and Cambon et al. (1981) for spectral closures. These closures are depending on a single numerical parameter only. Reflection of the difficulty of the task, many highly refined proposals for modeling this term have been proposed up to the intractable (see for instance [START_REF] Schiestel | Modeling and simulation of turbulent flows[END_REF]. Moreover, a recent study conducted by [START_REF] Mishra | Towards approximating non-local dynamics in single-point pressure strain correlation closures[END_REF] contributes to validate a closure approximation with a single numerical parameter similar to the previously cited models. In the case of subfilter quantities, considered here, the influence of the cutoff on this linear term is implicitly contained in the (P ij ) sf s term. The diffusion terms J ij is modeled assuming the well-known gradient law

J ij = ∂ ∂x m ν ∂(τ ij ) sf s ∂x m + c k k sf s ǫ (τ ml ) sf s ∂(τ ij ) sf s ∂x l (11)
where c k is a numerical coefficient.

∂ǫ ∂t + ∂ ∂x k (ū k ǫ) = c ǫ 1sf s ǫ k sf s P sf s -c ǫ 2sf s ǫ 2 k sf s + J ǫ (12)
The coefficient appearing in the destruction term of Eq. ( 12) is then given by (Chaouat and Schiestel, 2009; Chaouat, 2017)

c ǫ 2sf s = c ǫ 1 + ∆c ǫ [1 + βϑ α c ] γ (13)
where αγ = 2/3 and

β = [2/(3C K )] γ , C K is the Kolmogorov constant close to 1.5, ϑ c = κ c L where L
denotes the turbulence length-scale while κ c is the cutoff wave number computed using the grid size of the mesh ∆ as

κ c = π/∆, ∆c ǫ = c ǫ 2 -c ǫ 1 , c ǫ
c ǫ 1sf s = c ǫ 1 .
The term J ǫ denotes the diffusion term that is modeled using a tensorial gradient law.

J ǫ = ∂ ∂x j ν ∂ǫ ∂x j + c ǫ k sf s ǫ (τ jm ) sf s ∂ǫ ∂x m (14)
where c ǫ is a numerical coefficient. The subfilter model relying on Eqs. ( 6) and ( 12) has been also extended to low Reynolds numbers for approaching walls (see Chaouat and Schiestel, 2013).

Variance of the turbulent passive scalar field in spectral space

The new PITM development of the dissipation rate equation for the passive scalar variance is founded in spectral space in a way similar to the development of the dissipation rate equation for kinetic energy introduced in earlier works, and for this reason it will be considered first. The transported scalar is considered as passive with no direct action on the flow dynamics. It can represent temperature or any other transported contaminant and will be denoted θ. The main quantity representing the properties of turbulent scalar distribution is the variance of the scalar which is the counterpart of the kinetic energy of turbulence. The half variance k θ is defined as

k θ (X, t) = θ ′ θ ′ (X,
E θ (X, κ, t) = θ ′ θ ′ (X) ∆ (κ, t)/2
for the mean kinetic energy spectrum in homogeneous anisotropic turbulence and has been established in A according to previous works [START_REF] Tennekes | A first course of turbulence[END_REF][START_REF] Monin | Statistical Fluid Mechanics[END_REF]Schiestel, 1983aSchiestel, , 1983b;;[START_REF] Launder | Simulation of turbulence models and their applications[END_REF]Launder et al., 84b). This equation reads

∂E θ (X, κ, t) ∂t + u k (X) ∂E θ (X, κ, t) ∂X k = P θ (X, κ, t) + T θ (X, κ, t) +J θ (X, κ, t) -E θ (X, κ, t) (15)
where in the right hand side of this equation, P θ is the production of variance of the scalar by mean gradients of the scalar, T θ is the spectral transfer entrained by the eddying motions in the inertial cascade, J θ is the diffusion term and E θ denotes the dissipation term of the scalar variance. Without loosing generality, but for sake of simplification, we restrict provisionally the analysis to homogeneous flows so that the diffusion term vanishes and the variable X is omitted in the following. However, for practical non-homogeneous turbulent flows applications, the tangent homogeneous space approximation will be used. The production term P θ detailed in equation (74) of A involves the turbulent scalar fluxes ϕ iθ (X, κ, t) which have to be modeled carefully because they control the transport properties of turbulence. A physical approximation inspired from the GGDH hypothesis of Launder (1988) leads to

ϕ iθ (X, κ, t) = -ϕ ij (X, κ, t)τ (X, t) ∂ θ ∂x j (16)
where τ (X, t) is a characteristic time scale of turbulence. The hypothesis will be detailed further in the next section devoted to the physical space. The PITM method then uses a spectral decomposition where the wave number splitting can be viewed as a spectral filter [START_REF] Yoshizawa | A statistically-derived subgrid model for the large-eddy simulation of turbulence[END_REF]. Like in multiple-scale models (Schiestel, 1988), the spectral equation ( 15) can be partially integrated over the wave number range such as [0, κ c ], [κ c , κ e ] and [κ e , ∞[ where κ e denotes a wave number that is larger than κ c but it can be different from κ d . As a result, one can obtain the following approximate equations

∂k θles ∂t = P θ[0,κc] -F θ (κ c , t) (17)
∂k θsf s ∂t = P θ[κc,κe] -F θ (κ e , t) + F θ (κ c , t) (18)
0 = F θ (κ e , t) -ǫ θ[κe,∞[ (19)
where

k θsf s = κe κc E θ (κ, t)dκ , k θles = κc 0 E θ (κ, t)dκ , k θ = k θles + k θsf s (20)
F θ (κ c , t) = F θ (κ c , t) -E θ (κ c , t) ∂κ c ∂t (21)
F θ (κ e , t) = F θ (κ e , t) -E θ (κ e , t) ∂κ e ∂t (22)
and

ǫ θ[κe,∞] = ∞ κe E θ (κ, t)dκ (23)
where ǫ θ[κe,∞] ≈ ǫ θ . So that Eq. ( 18) can be rewritten as

∂k θ sf s ∂t = P θ[κc,κe] + F θ (κ c , t) -ǫ θ (24)
Equations ( 21) and ( 22) result from the derivative of the integral with variable bounds. The total variance transfer F θ (κ e , t) through the variable cutoff κ e is equal to the sum of the local spectral flux F θ (κ e , t) and the transfer due to the variation in the splitting wavenumber. At the wavenumber κ e , all the preceding hypotheses imply F θ (κ e ) ≈ ǫ θ , the turbulence Reynolds number being supposed to be large. For the scalar variance, the order of magnitude of the leading terms in Eq. ( 24) is

O(θ 2 u/l θ )
where l θ is the turbulence macroscale involved in the scalar physical processes. Consequently,

κ e -κ c = O(1/l θ ) = O(ǫ θ /θ 2 u) leading to the equation κ e -κ c = ζ θ ǫ θ k θsf s k 1/2 sf s (25)
where ζ θ is an adjustable coefficient chosen such that the spectral contribution of the variance beyond κ e is negligible. Equation ( 25) is the scalar counterpart of the dynamical equation already developed in the dynamical model [START_REF] Schiestel | Towards a new partially integrated transport model for coarse grid and unsteady turbulent flow simulations[END_REF][START_REF] Chaouat | A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows[END_REF]. Combining these equations together yields

∂ǫ θ ∂t = ǫ θ k θsf s ∂k θsf s ∂t + ǫ θ 2k sf s ∂k sf s ∂t + ǫ θ κ e -κ c F θ (κ e , t) -F θ (κ e , t) E θ (κ e , t) - F θ (κ c , t) -F θ (κ c , t) E θ (κ c , t) (26)
Using the transport equations for k sf s and ǫ recalled in section 2.4 and Eq. ( 24), one can obtain the resulting equation for the dissipation-rate ǫ θ written in a more compact form as

∂ǫ θ ∂t = c ǫ θθ 1 sf s P θsf s ǫ θ k θsf s + c ǫ θk 1 sf s P sf s ǫ θ k sf s -c ǫ θk 2 sf s ǫ θ ǫ k sf s -c ǫ θθ 2 sf s ǫ 2 θ k θsf s (27)
where

P θsf s = P θ[κc,κe] + F θ (κ c ) ( 28
)
c ǫ θθ 1 sf s = 1 , c ǫ θk 1 sf s = 1 2 , c ǫ θk 2 sf s = 1 2 (29)
and

c ǫ θθ 2 sf s = 1 - k θsf s (κ e -κ c )E θ (κ e ) F θ (κ e ) ǫ θ - F θ (κ e ) ǫ θ - E θ (κ e ) E θ (κ c ) F θ (κ c ) ǫ θ F θ (κ c ) F θ (κ c ) -1 (30)
Setting κ c ≪ κ e , E(κ d ) ≪ E(κ c ), and E θ (κ e ) ≪ E θ (κ c ), and also considering that F θ (κ e ) = ǫ θ as indicated by Eq. ( 19), Eq. ( 30) reduces then to

c ǫ θθ 2 sf s = 1 - k θsf s κ e E θ (κ e ) F θ (κ e ) ǫ θ -1 (31)
When κ c goes to zero, that is to say when the filter width in physical space goes to infinity in an homogeneous turbulence field (or locally homogeneous), one recovers the equation used in statistical RANS closure. Hence, the equation can be written as

∂ǫ θ ∂t = c ǫ θθ 1 P θ ǫ θ k θ + c ǫ θk 1 P ǫ θ k -c ǫ θk 2 ǫ θ ǫ k -c ǫ θθ 2 ǫ 2 θ k θ ( 32
)
where

c ǫ θθ 1 = 1 , c ǫ θk 1 = 1 2 , c ǫ θk 2 = 1 2 (33)
and

c ǫ θθ 2 = 1 - k θ κ e E θ (κ e ) F θ (κ e ) ǫ θ -1 (34)
Several modeling papers from the 1970's have proposed modeled forms of the dissipation-rate of the scalar. One of the most complete forms can be found in the work of [START_REF] Yoshizawa | Statistical modelling of passive-scalar diffusion in turbulent shear flows[END_REF] and also several authors such as [START_REF] Nagano | A two-equation model for heat transport in wall turbulent shear flows[END_REF] as well as [START_REF] Shikazono | Second-moment closure for turbulent scalar transport at various Prandtl numbers[END_REF]. In the earlier formulations of Eq. ( 32), it happened that some terms were obviously missing. But the form of Eq.

(32) obtained in the present case by a theoretical formalism seems indeed to be one of the most general considering the different source terms involved in this equation [START_REF] Newmann | Modelling the behaviour of homogeneous scalar turbulence[END_REF][START_REF] Jones | Closure of the Reynolds stress and scalar flux equations[END_REF]. Physically, the value of ǫ θ which can be viewed as the flux of variance in the inertial cascade must remain the same whatever the location of the spectral splitting. Comparing Eq. ( 31) with Eq. ( 34) leads to

c ǫ θθ 2 sf s = c ǫ θθ 1 + k θsf s k θ c ǫ θθ 2 -c ǫ θθ 1 ( 35
)
where the numerical value c ǫ θθ 1 supposed to take on the value unity can be slightly changed in practice.

As a result of interest, one can remark that Eqs. ( 27) and ( 32) have exactly the same general form as the usual modeled equation of the dissipation of variance found in the open literature and especially in Refs. [START_REF] Yoshizawa | Statistical modelling of passive-scalar diffusion in turbulent shear flows[END_REF][START_REF] Nagano | A two-equation model for heat transport in wall turbulent shear flows[END_REF] but also in Refs. [START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF]Mathieu et al., 1984b;Launder et al., 1984b;[START_REF] Ruffin | Investigation of characteristic scales in variable density turbulent jets using a second-order model[END_REF][START_REF] Gharbi | Numerical modelling of variable density turbulent jets[END_REF]. The development of the PITM method leads to find that it is the c ǫ θθ 2 sf s coefficient that controls the partitioning of scalar variance like the c ǫ 2sf s coefficient in the epsilon equation was controlling the kinetic energy. It is not an initial choice but the result of the development of the method.

4

Parametrization of the variable coefficients in the passive scalar dissipation equations

Molecular Prandtl numbers near unity

The ratio k θsf s /k θ appearing in Eq. ( 35) must be calibrated as a function of the location of the cutoff wavenumber. As indicated in Refs. [START_REF] Corrsin | Further generalizations of Onsager's cascade model for turbulent spectra[END_REF][START_REF] Tennekes | A first course of turbulence[END_REF] for Prandtl number near to unity, the spectrum of the scalar in the equilibrium range can be approximated by

E θ (κ) = C θ ǫ θ ǫ -1/3 κ -5/3 (36)
where C θ is a constant coefficient close to 0.5. The value of this numerical constant remains however approximate. Even to these present days, it is still not accurately well calibrated and is subject to discussion in the literature as pointed out by [START_REF] Favre | La turbulence en mécanique des fluides[END_REF]. The spectrum of the scalar θ given by Eq. ( 36) can be extended in the whole range domain of the wavenumbers using the very likely hypothesis

E θ (κ) = C θ ǫ θ C K ǫ E(κ) (37)
where E(κ) is the energy density spectrum of turbulence. Equation (37) must verify the limiting

condition lim κ→∞ E θ (κ) = C θ ǫ θ ǫ -1/3 κ -5/3
. The analytical integration made in B.1 yields the practical result

c ǫ θθ 2 sf s = c ǫ θθ 1 + ∆c ǫ θθ [1 + βϑ α c ] γ (38)
where ∆c ǫ θθ = c ǫ θθ 2 -c ǫ θθ 1 . As expected, the result is analogous to the formula previously obtained for the dynamical equations for kinetic energy. We mention also that for P r = 1, Eq. ( 37) reflects an analogy between the velocities and scalar fields corresponding to the well known Reynolds analogy.

Small molecular Prandtl numbers

This situation corresponds to the case of liquid metals that can have specific industrial applications. In this case, the inertial subrange of the variance spectrum is shorter and terminated by a steep inertial diffusive subrange. In this subrange, the spectral flux of variance of the scalar decreases drastically because of the dissipation of variance due to high molecular diffusivity. As indicated in Refs. [START_REF] Batchelor | Small scale variation of convected quantities like temperature in a turbulent fluid. Part 1. General discussion and the case of small conductivity[END_REF]; Tennekes and Lumley, 1972) the spectrum of the scalar variance can be represented by the function

E θ (κ) = C θ ǫ θ ǫ -1/3 κ -5/3 exp - 3 2 C θ (κη θ ) 4/3 (39)
with the scalar microscale defined by η θ = (σ 3 /ǫ) 1/4 where C θ = 1.5 is a constant coefficient such that r . The spectrum shape associated with Eq. ( 39) is shown in Fig. 3. Dealing with an exponential factor results in cumbersome formulas in the model equations. So, in practice a simpler approach is retained. Equation ( 39) is then replaced by

C θ = k θ C K ǫ/(kǫ θ )
LogE k k k c e k k SGS LES F c F e O e q q q q inertial convective subrange (-5/3) inertial diffusive subrange production zone k d ( ) PrRe L 3/4 t k = LogE LogE q H q q k H
E θ (κ) = C θ ǫ θ ǫ -1/3 κ -5/3 H(κ H -κ) (40)
where κ H = 1/η θ , and H is the Heaviside function implying that E θ (κ) = 0 for κ ≥ κ H . The spectral vanishing value of wavenumber is then obtained for κη θ = 1. So that the dimensionless variable ϑ is dropping for ϑ H = (P r Re t ) 3/4 where Re t = k 2 /νǫ denotes the turbulent Reynolds number. This dropping value can be expressed equivalently as ϑ H = k 3/2 /(η θ ǫ) = (σk 2 /ǫ) 3/4 . Physically, the dimensionless group P r Re t is interpreted like the turbulent Peclet number denoted P e t = P r Re t . Obviously, as seen in Fig. 4, the value of κ e can be larger or smaller than η -1 θ . The spectrum [START_REF] Mishra | Towards approximating non-local dynamics in single-point pressure strain correlation closures[END_REF] valid in the inertial law can be extended in the whole wavenumber range by writing

E θ (κ) = C θ ǫ θ C K ǫ E(κ)H(κ H -κ) (41)
where E(κ) is the spectrum valid in the whole wave number range given by Eq. ( 83). The exact final expression of the coefficient c ǫ θθ 2 sf s is developed in B.2 and reads

c ǫ θθ 2 sf s =            c ǫ θθ 1 + ∆c ǫ θθ [1 + βϑ α c ] -γ -1 + βP e 3α/4 t -γ 1 -1 + βP e 3α/4 t -γ (ϑ c < ϑ H ) c ǫ θθ 1 (ϑ c > ϑ H ) (42)

Large molecular Prandtl numbers

This situation corresponds to the case of poorly conducting fluids or high viscous fluids like most of oils.

In this case, the inertial subrange is followed by a viscous-convective subrange with a negative slope of minus unity and finally a viscous-diffusive subrange in which the spectrum undergoes strong decay as indicated in Refs. [START_REF] Batchelor | Small scale variation of convected quantities like temperature in a turbulent fluid. Part 1. General discussion and the case of small conductivity[END_REF][START_REF] Tennekes | A first course of turbulence[END_REF]). In the viscous-convective subrange, the scalar fluctuations are broken and decreased by the strain-rate field but the scalar diffusivity is not yet active. The shape of the spectrum has been studied by [START_REF] Batchelor | Small scale variation of convected quantities like temperature in a turbulent fluid. Part 1. General discussion and the case of small conductivity[END_REF]. For the wave number κ ≥ 1/η K , it can be shown that the viscous convective subrange of the spectrum is of the form

E θ (κ) = c θ ǫ θ ν ǫ 1/2 κ -1 (43)
where c θ is a constant coefficient. The viscous convective subrange is followed by the viscous-diffusive subrange which is characterized by the role of scalar diffusivity acting on very small scales. In this region, the spectrum takes on the form

E θ (κ) = c θ ǫ θ ν ǫ 1/2 κ -1 exp -c θ (κη * θ ) 2 (44)
where

η * θ = η K (σ/ν) 1/2 = η K / √ P r
is the smallest scale of the viscous-diffusive subrange and c θ is a constant coefficient. The corresponding wave numbers are then computed as κ K = 1/η K and κ S = 1/η * θ . The junctions between the different subranges occur for κ = κ K and κ = κ S respectively preserving continuity. In particular, for κ = κ K , the spectrum E θ (κ K ) given by Eq. ( 36) of the scalar in the equilibrium range with a slope κ -5/3 is equal to the spectrum E θ (κ K ) given by Eq. ( 43) of the viscousconvective subrange with a slope κ -1 so that it is found that c θ = C θ ≈ 1.5. The dimensionless wave numbers ϑ = κL associated with the Kolmogorov scale η K and the smallest scale η * θ are ϑ K = Re by Eq. ( 44) is replaced by a simple form as

E θ (κ) = c θ ǫ θ ν ǫ 1/2 κ -1 H(κ S -κ) (45)
implying that E(κ) = 0 for κ ≥ κ S . In the wave number range [0, κ K ], the spectrum still evolves according to Eq. ( 83). The wave number range [0, κ S ] is then decomposed into two wave number ranges introducing the cutoff wave number κ c where κ c < κ S or κ c > κ S . In the first wave number range [0, κ K ], the spectrum E θ 1 (κ) is defined as

E θ 1 (κ) = ξ k θ k E(κ) (46)
where E(κ) is given by Eq. ( 83) whereas in the second domain [κ K , ∞[, the spectrum E θ 2 (κ) is deduced from Eq. ( 45)

E θ 2 (κ) = ξc θ ǫ θ ν ǫ 1/2 κ -1 H(κ S -κ) ( 47
)
where ξ is a coefficient of normalization. After an analytical integration made in B.3, we obtain the exact expression for the c ǫ θθ 2 sf s coefficient

c ǫ θθ 2 sf s =                          c ǫ θθ 1 + ∆c ǫ θθ [1 + βϑ α c ] -γ -1 + βRe 3α/4 t -γ + c θ ǫ θ k θ ν ǫ 1/2 ln P 1/2 r 1 -1 + βRe 3α/4 t -γ + c θ ǫ θ k θ ν ǫ 1/2 ln P 1/2 r (ϑ C < ϑ K ) c ǫ θθ 1 + ∆c ǫ θθ c θ ǫ θ k θ ν ǫ 1/2 ln ϑ S ϑc 1 -1 + βRe 3α/4 t -γ + c θ ǫ θ k θ ν ǫ 1/2 ln P 1/2 r (ϑ K < ϑ C < ϑ S ) c ǫ θθ 1 (ϑ S < ϑ C ) (48)
Introducing the dimensionless ratio R that characterizes the turbulence scale of the passive scalar field compared to the dynamic scales [START_REF] Kenjeres | Contribution to elliptic relaxation modelling of turbulent natural and mixed convection[END_REF][START_REF] Kenjeres | Capturing transient effects in turbulent flows over complex urban areas with passive pollutants[END_REF]) R = (k θsf s ǫ)/(k sf s ǫ θ ) then, the factor 96) can be approximated in a simple way using the relations

ζ = (ε θ /k θ )(ν/ε) 1/2 involved in equation (
(ε θ /k θ ) = ε/(kR) and (ν/ε) 1/2 = k/ εRe 1/2 t
, thus implying

ε θ k θ ν ε 1/2 = 1 R Re -1/2 t ( 49
)
At a first sight, it appears that these analytical coefficients seem fairly complicated. But in fact, they can be easily implemented in CFD codes, because the extra computing time is entirely negligible. As emphasized above, this is the price to pay for dealing with an extended Prandtl number range.

5

Final equation set for passive scalar turbulence

Turbulent passive scalar field equations in physical space

The transport equation of the subfilter scale scalar variance including additional turbulent diffusion terms in non-homogeneous turbulence analogous to Eq. ( 24) reads

∂k θ sf s ∂t + ∂ ∂x k ūk k θsf s = P θsf s + J θ -ǫ θ (50)
The production term P θsf s is given by

P θsf s = -(τ jθ ) sf s ∂ θ ∂x j (51)
where (τ iθ ) sf s represents the turbulent scalar flux vector which is of primary importance for the determination of turbulent transfer properties. In a second order closure framework, it is advisable to solve the transport equation for (τ iθ ) sf s but all the problems mentioned previously for the linear part of the pressure-strain correlation are still posed here for the pressure-scalar gradient correlation term (Π iθ ) sf s [START_REF] Launder | Simulation of turbulence models and their applications[END_REF][START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF][START_REF] Gatski | Second-moment and scalar flux representations in engineering and geophysical flows[END_REF]. However, in the present case, as mentioned above, due focus is given to the variance of scalar fluctuation and it is assumed in a preliminary step

(τ iθ ) sf s = -c τ θ (τ im ) sf s k sf s ǫ ∂ θ ∂x m (52)
where c τ θ is a numerical coefficient. This modeling is in the present case here also strengthened by the fact that the turbulent stress tensor τ ij appearing in Eq. ( 52) is computed by a second-moment closure (SMC) which allows to accurately reproduce the turbulence anisotropy and not by simple eddy viscosity models. So, the tensorial diffusivity benefits obviously from the second order closure in turbulent stresses. This is an important point because the turbulent scalar fluxes stand as input in the scalar variance equation and also its dissipation rate. Comparisons of turbulent scalar fluxes to known DNS calculations are given in the following in order to ensure that proper scalar flux predictions are obtained using Eq. ( 52). As the turbulent stresses are computed in the framework of second moment closure, the diffusion term J θ is modeled also using a well-known tensorial gradient law hypothesis

J θ = ∂ ∂x m σ ∂k θ sf s ∂x m + c k θ (τ mj ) sf s k sf s ǫ ∂k θ sf s ∂x j (53)
where c k θ is a constant coefficient. The final equation of the scalar variance dissipation-rate ǫ θ is deduced from Eq. ( 27) including the additional diffusion term as follows

∂ǫ θ ∂t = c ǫ θθ 1 sf s P θsf s ǫ θ k θsf s + c ǫ θk 1 sf s P sf s ǫ θ k sf s -c ǫ θk 2 sf s ǫ θ ǫ k sf s -c ǫ θθ 2 sf s ǫ 2 θ k θsf s + J ǫ θ ( 54
)
where the diffusion term J ǫ θ is modeled assuming the tensorial gradient law hypothesis

J ǫ θ = ∂ ∂x m σ ∂ǫ θ ∂x m + c ǫ θ (τ mj ) sf s k sf s ǫ ∂ǫ θ ∂x j (55)
and c ǫ θ is a constant coefficient. The coefficients c ǫ θθ 1 sf s , c ǫ θk 1 sf s , c ǫ θk 2 sf s and c ǫ θθ 2 sf s appearing in Eq. ( 54) are given by Eqs. ( 29) and Eq. ( 84) for P r ≈ 1, Eq. ( 89) for P r ≪ 1 and Eq. ( 96) for P r ≫ 1.

Practical choice of the spectrum parameters and model coefficients

As it was mentioned in section 4, the energy density spectrum is governed by the parameters α and γ.

Trying to stick to the values suggested by the PITM study, we use here the values

c ǫ θθ 1 = 1, c ǫ θk 1 = 1/2, c ǫ θk 2 = 1/2,
c ǫ 1 c ǫ 2 c k c ǫ c ǫ θθ 1 c ǫ θk 1 c ǫ θk 2 c ǫ θθ 2 c τ θ c k θ c ǫ θ 1.

6

Application to the fully developed turbulent channel flow subjected to a constant heat flux

Filtered equations for the simulation in physical space

The filtering of the instantaneous Navier-Stokes equations gives rise to the mass conservation and resolved scale motion equations as ∂ ūj ∂x j = 0 (56)

∂ ūi ∂t + ∂(ū i ūj ) ∂x j = - 1 ρ ∂ p ∂x i + ν ∂ 2 ūi ∂x j ∂x j - ∂(τ ij ) sf s ∂x j (57)
where (τ ij ) sf s = u i u j -ūi ūj . In Eqs. ( 56) and ( 57), the commutation terms due to the fact that in the general case the filtering operation does not commute with the space or time derivative have been neglected (Chaouat and Schiestel, 2013;Chaouat, 2017c). The filtered transport equation for the passive scalar is

∂ θ ∂t + ∂( θū j ) ∂x j = σ ∂ 2 θ ∂x j ∂x j - ∂(τ jθ ) sf s ∂x j (58)
where (τ iθ ) sf s = u i θ -ūi θ. This set of equation is solved to carry on the simulation. It is coupled with the PITM model to account for subfilter turbulence.

Plane channel flow subjected to a constant heat flux on both walls

This section aims to illustrate the theoretical developments accounting for passive scalar transport through basic numerical simulations. We perform the fully developed turbulent channel flow subjected to uniform heat fluxes q w at both walls as indicated in Fig. 6 for the Reynolds number R τ = u τ δ/ν = 395 where u τ denotes the friction velocity and δ is the half channel width, and for three values of the molecular Prandtl number P r = 0.1, 1 and 10. Physically, these Prandtl numbers are chosen as representative of heat transfer in liquid metals such as mercury or liquid sodium P r ≈ 0.015), gases (P r ≈ 1), water (P r ≈ 5 -7) and seawater (P r ≈ 13 at 0 • C). We have carried out a RANS-LES simulation based on PITM formalism of the turbulent plane channel flow in order to validate the subfilter model equations including both the dynamic equations and the thermal equations. In order to focus into the passive scalar field that constitutes the original part of this present work, the purely dynamical illustrative results are not presented here but they can be found in previous papers of the authors themselves (Chaouat andSchiestel 2005, 2009). The PITM simulations are then compared with reference DNS performed by Chaouat (2018), Chaouat and Peyret (2019).

Numerical procedure

The simulations are performed on several meshes of coarse and medium grid resolutions to study the grid effect on the scalar variable, turbulent heat fluxes (τ iθ ) sf s , and variance of the passive scalar k θ with emphasis placed on the sharing out between the modeled and resolved scales. The dimensions of the channel in the streamwise, spanwise and normal directions along the axes and Prandtl number values studied here, the grids resolution in the (x 1 , x 2 , x 3 ) directions as well and the grid spacings in wall unit are listed in Table 3. The meshes are uniform in the streamwise and spanwise directions (x 1 , x 2 ) and the grid spacings in wall unit are ∆ + 1 = ∆ + 2 = 60 for the coarse mesh and ∆ + 1 = ∆ + 2 = 30 for the medium mesh, respectively. The spacing in the center of the channel along the normal direction to the wall is ∆ + 3c = 16. In this direction x 3 , the grid points are distributed in different spacings with a refinement near the wall to compute accurately both the velocity and thermal boundary layers according to the transformation x 3j = tanh [ξ j F (ξ j ) atanh a]/2 where ξ j = -1+2(j-1)/(N 3 -1)

x 1 , x 2 , x 3 are respectively O 1 X 2 X 3 U 1 X q w q w δ 2
(j = 1, 2, • • • N 3 ), F (ξ j ) = (1 + ξ 2 j )/2
, the parameter a is a coefficient set to 0.990 for N 3 = 84 and 0.9875 for the N 3 = 128 leading to ∆ + 3w = 0.2. For both cases, it is found that the larger grid spacing in the normal direction is then smaller than 16 in wall unit ∆ + 3 < 16. This choice of grid size allows generally to get a significant part of the scalar variance in the subfilter part. However in the case of small Prandtl numbers, the vanishing wavenumber of the spectrum κ H is greatly diminished, inescapable clue of high scalar dissipation implying a very smaller subfilter part. But this is a physical consequence and does not mean at all that the grid is too fine. In terms of an order of resolution purely for information purposes, the grid-point and the spacings associated with the direct numerical simulation performed by [START_REF] Chaouat | Investigation of the wall scalar fluctuations effect on passive scalar turbulent fields at several Prandtl numbers by means of direct numerical simulations[END_REF] are also indicated in Table 3. As a result of interest, the ratio of the DNS grid-points over the PITM grid-points R(P r ) = N DN S /N P IT M varies from 113 for P r = 0.1 to 2378 for P r = 10. These values clearly show that all PITM grids are of coarse resolution in comparison with DNS grids and also highly resolved LES grids. It can be emphasized that the resolution of the grid for a direct numerical simulation of scalar turbulence is loose at low Prandtl number but fine at large Prandtl number according to the power law of the Prandtl number N ∝ R that the present grids used here, especially at the Prandtl number P r = 10 are in fact very coarsened in comparison with the DNS grids for such PITM simulations that run in LES mode far away the walls because much more scales need to be modeled. The simulations are performed using the numerical code developed by [START_REF] Chaouat | An efficient numerical method for RANS/LES turbulent simulations using subfilter scale stress transport equations[END_REF] which is based on the finite volume technique and optimized with message passing interface (MPI). The equations are integrated in time by means of an implicit Runge-Kutta scheme of fourth-order accuracy in time. The Navier-Stokes equations are solved in space by a centered scheme of fourth-order accuracy in space whereas the transport equations of turbulence and scalar fields are solved by an upwind scheme of second-order accuracy in space. Equations ( 58) and ( 52) associated with the mean scalar variable are independent from the set of Eqs. ( 50) and ( 54) associated with the scalar variance. A constant pressure gradient term is included in the momentum equation to balance the viscous effects at the walls and an extra source term is added in the transport equation of the scalar variable aiming to get a periodic condition between the inlet and outlet plane sections of the channel [START_REF] Chaouat | DNS of passive scalar transport fields in turbulent flow at low and high Prandtl numbers[END_REF][START_REF] Chaouat | Investigation of the wall scalar fluctuations effect on passive scalar turbulent fields at several Prandtl numbers by means of direct numerical simulations[END_REF]. The boundary conditions imposed at the lower and upper walls located at x 3 = 0 and 2δ of the channel are no slip velocity conditions ūi = 0. A constant heat flux corresponding to an isoflux boundary condition given by q w = -κ(∂ θ/∂x 3 ) w is applied at the upper and lower walls.

Numerical results

As usually, the results are presented in dimensionless form using the half channel width δ, the wall friction velocity u τ and the friction scalar variable θ w = q w /(ρc p u τ ) where c p denotes the specific heat at constant pressure. The transformed variable Θ + = θ + w -θ + is considered in order to analyze the present results because the mean statistical quantities including the velocity and normal turbulent stresses are symmetric about the centerline of the channel. Figure 7 describes the profile of the mean scalar variable Θ + versus the logarithmic wall distance for several Prandtl numbers returned by the two simulations performed on the coarse and medium meshes, where the brackets . denotes the averaging in homogeneous directions of the flow. As known, the logarithmic region gets more extended when the Prandtl number increases from 0.1 to 10. An excellent agreement is observed between the DNS and PITM simulations at each Prandtl number confirming that the thermal boundary layer is accurately resolved in the wall region. As a preliminary requirement, it is necessary to verify the predictions of the turbulent heat fluxes. Indeed these turbulent fluxes not only control the mean distribution of temperature but also appear in the driving source term of the variance, they are thus twice important.

Figures 8 and9 exhibit the profiles of the streamwise and normal turbulent heat fluxes |q 1 | = u ′+ 1 θ ′+ and q 3 = u ′+ 3 θ ′+ , respectively, versus the wall distance for the Prandtl numbers P r = 0.1, 1 and 10 and for the two simulations performed on the coarse and medium meshes. Overall, it is found that the PITM simulation provides turbulent heat fluxes in good agreement with the reference DNS data. The intensity of the streamwise turbulent heat flux for the Prandtl number P r = 0.1 and unity appears slightly overpredicted when moving towards the center of the channel because of the resolved contribution which decreases a little bit too slowly away from the walls. This is probably due to the coarse grid resolution in this flow region. As expected, the normal turbulent heat flux is accurately predicted for each Prandtl number according to the DNS data. These figures show that the subfilter turbulent heat fluxes are relatively high in the near wall region revealing the presence of a thermal peak for P r = 1 and 10 followed afterward by a rapid decrease when going away from the wall, while the resolved turbulent heat fluxes increase and decrease more slowly versus the wall distance. The curves associated with the resolved turbulent heat fluxes present a regular evolution from the wall to the centerline of the channel. Because of the grid resolution that has been chosen, the subfilter part of the turbulent heat flux as determined by the location of the cutoff wave number κ c in the spectrum of the scalar variable E θ (κ) appears of lower intensity than its corresponding resolved part almost everywhere in the channel even if the sharing out of the energy is appreciably modified with respect to the Prandtl number. Indeed, as the Prandtl number increases from the lower value P r = 0.1 to the higher value 10, the subfilter part is increasing relatively to its resolved part and conversely, the resolved part is decreasing. This is an important point suggesting that the model is consistent between the two extremes limits of the energy spectra splittings of E(κ) and E θ (κ) that are RANS and DNS. As expected, it is found that the subfilter heat flux energy increases as the grid-size increases since larger scales must be modeled and vice versa, the resolved heat flux energy decreases since less scales are computed. But as a remarkable result, it is found that the sum of each contribution, modeled and resolved parts for the coarse and medium corresponds to the DNS fairly well. This observation appears so clear for the normal turbulent heat fluxes q 3 = u ′+ 3 θ ′+ but is less marked for the streamwise turbulent heat flux q 1 = u ′+ 1 θ ′+ . Hence, the analysis of the results with respect the two grids considered confirms that the PITM method returns heat fluxes in good agreement with DNS independently of the spectral cutoff location. It is worth mentioning that the mean scalar variable θ deduced from the instantaneous filtered variable θ solution of Eq. ( 58) is here governed by the normal subfilter turbulent heat (τ 3θ ) sf s .

As the mean scalar profile shown in Fig. 7 is fairly well predicted from the simulation at each Prandtl number, it is not at all surprising that the normal heat flux profiles agree very well with the DNS data for all grids. So, this outcomes demonstrates that the closure (52) used for the turbulent heat flux assuming the generalized gradient diffusion hypothesis GGDH [START_REF] Daly | Transport equations in turbulence[END_REF] is entirely satisfactory for solving the scalar-velocity correlations in the present application. This one is generally successful in heat transfer of practical engineering and geophysical flows when the Reynolds stress anisotropy is reasonably well captured both in RANS and LES [START_REF] Launder | On the computation of convective heat transfer in complex turbulent flows[END_REF][START_REF] Suga | Improvement of second-moment closure for turbulent obstacle flow and heat transfer[END_REF][START_REF] Kenjeres | LES,T-RANS and hybrid simulations of thermal convection at high Ra number[END_REF][START_REF] Kenjeres | Capturing transient effects in turbulent flows over complex urban areas with passive pollutants[END_REF]. As regards temperature variance predictions, figure 10 displays the root-mean square of the scalar variance θ + rms = θ ′+ θ ′+ = √ 2k θ computed by means of the subgrid and resolved fluctuations as follows

θ ′ θ ′ = [ θθ -θ θ ] + [ θ θ -θ θ ] = 2 k θsf s + k θles (59)
respectively, for the simulations performed on the coarse and medium meshes. A quantitative agreement is obtained with the DNS although the PITM is performed on a relatively coarse grid, especially in the streamwise and spanwise directions of the mesh. For Pr=0.1, the subfilter part of the scalar variance appears to be of lower intensity than the resolved part everywhere in the channel. This result is not a shortcoming but must be attributed to Eqs. (84), ( 89) and (96) accounting for the dimensionless numbers ϑ c , ϑ H and ϑ K . Indeed at very low Prandtl numbers, the variance spectrum is falling down very early and thus contains very few small eddies of the scalar. In other words, the high diffusivity efficiently cancel these small eddies. As has been observed previously for the turbulent heat fluxes, when the grid-size increases in space ∂∆(x i )/∂x i > 0, then a part of energy k θ contained into the resolved scales is removed and fed into the modeled spectral zone, whereas on the contrary, when ∂∆(x i )/∂x i < 0, a part of energy coming from the modeled zone is injected into the resolved scales but the total variance k θ remains almost the same for the two grids considered as it was also obtained for the turbulence energy k (Chaouat and Schiestel 2005, 2012; Chaouat, 2017c). Figure 11 shows the contours plots of the instantaneous filtered scalar field Θ for the Prandtl number P r = 1 in the (x 1 , x 3 ) mid-plane of the channel. The unsteady character of the scalar field is well visible. In particular, this figure puts in light the detachment of vortex in the normal direction to the wall according to the direct numerical simulation performed by Chaouat and Peyret, 2019 (see Figure 14) even if some details are lost in the present case because of the very coarse grid resolution. In this context, it can be recalled that RANS computations cannot reproduce instantaneous structures in time and space because of the RANS physical foundations relying in practice on a long-time averaging which is large in comparison with the turbulence time scale (Chaouat, 2017). As known, RANS determine only mean organized structures as for instance the parietal vortex shedding that develop in a channel flow with wall mass injection [START_REF] Chaouat | Reynolds stress transport modelling for steady and unsteady channel flows with wall injection[END_REF].

Table 3: Simulation parameters for the fully developed turbulent channel flow, grid-points N i , total number of grid-points N, grid resolution ∆ + i , in the streamwise, spanwise and wall-normal directions, x 1 , x 2 , x 3 , respectively, and Batchelor length-scale η θ relatively to the Kolmogorov length-scale η κ for the Prandtl numbers P r = 0.1, 1 and 10. Domain size L 1 × L 2 × L 3 = 6.4δ × 3.2δ × 2δ. Present PITM simulations performed at the Reynolds number R τ = u τ δ/ν = 395. DNS [START_REF] Chaouat | Investigation of the wall scalar fluctuations effect on passive scalar turbulent fields at several Prandtl numbers by means of direct numerical simulations[END_REF].

Case N 1 N 2 N 3 N(10 6 ) ∆ + 1 ∆ + 2 ∆ + 3c η θ DNS (Pr=0.

Concluding remarks

The present work was devoted to the extension of the PITM method initially developed for simulating unsteady turbulent flows to the associated passive scalar transport fields from a physical and analytical standpoint. Emphasis has been put on the methodology. In this framework, we make use of the concept of tangent homogeneous space to derive the spectral transport equation for the variance of the passive scalar. Using partial integration in the spectral space and spectral splitting techniques, we have then derived the equations for the subfilter variance k θsf s and its dissipation-rate ǫ θ in the physical space. As a result, it appears that the dissipation-rate equation for ǫ θ takes the same form as the The instantaneous transport equation for a passive scalar can be written as

∂θ ∂t + u j ∂θ ∂x j = σ ∂ 2 θ ∂x j ∂x j (60)
where σ = λ/ρc p = ν/P r is the thermal diffusivity, λ is the thermal conductivity, c p is the specific heat at constant pressure, ρ is the density, P r is the Prandtl number and ν denotes the molecular viscosity.

Using the decomposition θ = θ + θ ′ , the transport equation for the fluctuating temperature θ ′ is then

∂θ ′ ∂t + u j ∂θ ′ ∂x j = -u ′ j ∂ θ ∂x j -u ′ j ∂θ ′ ∂x j + σ ∂ 2 θ ′ ∂x j ∂x j + ∂ ∂x j u ′ j θ ′ (61)
We consider the general case of nonisotropic inhomogeneous turbulence. In this case, the two-point

temperature correlation φ θθ = θ ′ A θ ′ B (x A , x B )
, where θ ′ A is treated as a function of x A whereas θ ′ B is a function of x B , considering that x A and x B are the primary independent variables. The transport equation for φ θθ is

∂φ θθ ∂t + u jA ∂φ θθ ∂x jA + u jB ∂φ θθ ∂x jB = -u ′ jA θ ′ B ∂ θ A ∂x jA -θ ′ A u ′ jB ∂ θ B ∂x jB - ∂ ∂x jA θ ′ A u ′ jA θ ′ B - ∂ ∂x jB θ ′ A u ′ jB θ ′ B +σ ∂ 2 φ θθ ∂x jA ∂x jA + ∂ 2 φ θθ ∂x jB ∂x jB (62)
New independent variables defined by the vector difference ξ = x Bx A and the midway position X = 1 2 (x A + x B ) are then introduced in the present derivation in order to distinguish the effects of distance separation from the effects of space location. So that each variable can be regarded as a new function of the two variables ξ and X. The formula for derivatives in the change of variables are detailed in Ref. [START_REF] Hinze | Turbulence[END_REF], the main ones being

∂ ∂x k A = 1 2 ∂ ∂X k - ∂ ∂ξ k (63) and ∂ ∂x k B = 1 2 ∂ ∂X k + ∂ ∂ξ k (64)
We introduce also the tensors

φ jθ = u ′ jA θ ′ B and φ θj = θ ′ A u ′ jB .
Taking into account these considerations, the complete dynamic equation (62) for the double temperature correlation for incompressible fluid flow is rewritten as follows

∂φ θθ (X, ξ, t) ∂t + 1 2 ( u jA + u jB ) (X, ξ, t) ∂φ θθ (X, ξ, t) ∂X j = -φ jθ (X, ξ, t) ∂ θ A ∂x jA = -φ θj (X, ξ, t) ∂ θ B ∂x jB -u jB -u jA (X, ξ, t) ∂φ θθ (X, ξ, t) ∂ξ j - 1 2 ∂ ∂X j θ ′ A u ′ jA θ ′ B + θ ′ A u ′ jB θ ′ B (X, ξ, t) - ∂ ∂ξ j θ ′ A u ′ jB θ ′ B -θ ′ A u ′ jA θ ′ B (X, ξ, t) + σ 2 ∂ 2 φ θθ ∂X j ∂X j (X, ξ, t) + 2σ ∂ 2 φ θθ ∂ξ j ∂ξ j (X, ξ, t) (65)
Before taking the Fourier transform of this equation, it is convenient to introduce Taylor series expansion in space of the terms involving the velocities leading to

( u kB -u kA )(X, ξ) = ξ m ∂ u k ∂X m (X) + O(ξ 2 ), (66)
( u kA + u kB )(X, ξ) = 2 u k (X) + ξ m ξ p 2 ∂ 2 u k ∂X m ∂X p (X) + O(ξ 3 ), (67)
It is also useful to develop the production term involving the derivative of the mean temperature as follows

φ jθ ∂ θ A ∂x jA + φ θj ∂ θ B ∂x jB = φ jθ ∂ θ ∂X j + φ θj ∂ θ ∂X j - ξ m 2 ∂ 2 θ ∂X m ∂X j (φ jθ -φ θj ) + O(ξ 2 ) (68)
Here, we restrict the development only to its first term to get a tractable equation. This means that the Fourier terms are identical to the ones in homogeneous anisotropic turbulence. Using Eqs. (66), (67) and (68), Eq. (65) including contributions of homogeneous and nonhomogeneous terms becomes

∂φ θθ (X, ξ, t) ∂t + u j (X) ∂φ θθ (X, ξ, t) ∂X j = -φ jθ (X, ξ, t) ∂ θ A ∂X j -φ θj (X, ξ, t) ∂ θ B ∂X j -ξ m ∂ u j ∂X m ∂φ θθ ∂ξ j (X, ξ, t) - 1 2 ∂ ∂X j θ ′ A u ′ jA θ ′ B + θ ′ A u ′ jB θ ′ B (X, ξ, t) - ∂ ∂ξ j θ ′ A u ′ jB θ ′ B -θ ′ A u ′ jA θ ′ B (X, ξ, t) + σ 2 ∂ 2 φ θθ ∂X j ∂X j (X, ξ, t) +2σ ∂ 2 φ θθ ∂ξ j ∂ξ j (X, ξ, t) (69)
The nonhomogeneous terms that appear in Eq. (69) correspond to the usual terms in one-point equation whereas the other terms involving the distance ξ can be treated as in homogeneous anisotropic turbulence. So, this method can be viewed as considering the tangent homogeneous anisotropic field at the point X of the nonhomogeneous field.

A.2 Transport equation of the Fourier transform of the two-point correlation tensor

Considering the Fourier transform of φ θθ (X, ξ, t) that is expressed as

φ θθ (X, κ, t) = ∞ -∞ ∞ -∞ ∞ -∞ φ θθ (X, ξ, t) exp (-jκξ) dξ, (70)
the Fourier transform of the transport equation of the double velocity correlation in locally tangent spectral space then reads

∂ φ θθ (X, κ, t) ∂t + u j (X) ∂ φ θθ (X, κ, t) ∂X j = -φ jθ (X, κ, t) ∂ θ A ∂X j -φ θj (X, κ, t) ∂ θ B ∂X j +κ j ∂ u j ∂X m ∂ φ θθ ∂κ m (X, κ, t) - 1 2 ∂ ∂X j (S θj,θ + S θ,jθ ) -jκ j (S θ,jθ -S θj,θ ) + σ 2 ∂ 2 φ θθ ∂X j ∂X j (X, κ, t) -2σκ 2 φ θθ (X, κ, t) (71)
where

S θj,θ = θ ′ A u ′ jA θ ′ B and S θ,jθ = θ ′ A u ′ jB θ ′ B .

A.3 Transport equation of the spherical mean of the Fourier transform of the two-point correlation tensor

Considering φ θθ (X, t) and its Fourier transform φ θθ (X, κ, t), we define the spherical mean of the Fourier transform by the relation

ϕ θθ (X, κ, t) = [φ θθ (X)] ∆ (κ, t) = 1 4πκ 2 π 0 2π 0 φ θθ (X, κ, t)κ 2 sinθdθdφ (72)
The transport equation of the one-dimensional spectral tensor of the double temperature correlations is obtained by taking the Fourier transform and mean integration over spherical shells of Eq. (71). As mentioned in Refs. [START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF][START_REF] Cambon | Spectral modelling of homogeneous non-isotropic turbulence[END_REF], spherical averages allow to make some useful simplifications in the spectral equations. The directional information is lost but the spectral averaged correlation of the fluctuating temperatures is then only function of the wavenumber and not anymore of the wavevector. For sake of clarity, we denote E θ = ϕ θθ /2. In this case, the mistress transport equation for E θ then reads

∂E θ ∂t (X, κ, t) + u k (X) ∂E θ (X, κ, t) ∂X k = P θ (X, κ, t) + T θ (X, κ, t) + J θ (X, κ, t) -E θ (X, κ, t) (73)
where on the right hand side of Eq. ( 73), P θj is the production term defined by

P θ (X, κ, t) = -ϕ jθ (X, κ, t) ∂ θ ∂X j ( 74
)
T θ is the transfer term including two different contributions of different meaning

T θ (X, κ, t) = - 1 2 ξ m ∂φ θθ ∂ξ j (X, ξ, t) ∆ ∂ u j ∂X m - 1 2 ∂ ∂ξ j (S θ,jθ -S θj,θ (X, ξ, t)) ∆ (75)
J θ denotes the diffusion term given by

J θ (X, κ, t) = - 1 4 ∂ ∂X j S ∆ θ,jθ + S ∆ θj,θ + σ ∂ 2 E θ ∂X j ∂X j (X, ξ, t) (76)
and finally, E θ is the dissipation rate including two different contributions as

E θ = σ 2 ∂ 2 E θ ∂X j ∂X j (X, ξ, t) + 2σκ 2 E θ (X, ξ, t) (77)
Indeed, from the derivatives formula, it is straightforward to get

∂ 2 φ θθ (x A , x B ) ∂x jA ∂x jB ∆ = 1 4
∂ 2 ϕ θθ ∂X j ∂X j (X, ξ) - ∂ 2 ϕ θθ ∂ξ j ∂ξ j (X, ξ) (78)
showing the direct physical interpretation of the scalar dissipation term. As a result of interest, it is simple to see that the full integration in the spectral space of E θ when ξ goes to zero allows to recover the usual dissipation-rate ǫ θ The scalar variance k θ is obtained by integrating the spectrum (37) from 0 to infinity leading to

ǫ θ = ∞ 0 σ 2 ∂ 2 E θθ ∂X j ∂X j + 2σκ
k θ = ∞ 0 C θ ǫ θ C K ǫ E(κ)dκ = C θ ǫ θ k C K ǫ (81)
So that the reduced spectrum E * θ (ϑ) defined by E * θ (ϑ) = E θ (κ)/k θ L takes the simple expression written as E * θ (ϑ) = E(κ)/kL or equivalently E * θ (ϑ) = E * (ϑ). One can then derive the practical formulation of the dissipation equation of scalar variance for molecular Prandtl numbers near unity corresponding to

where ∆c ǫ θθ = c ǫ θθ 2 -c ǫ θθ 1 . As expected, the spectral properties of the scalar variance at Prandtl numbers near unity is very similar to the ones for velocity variance, expressing some Reynolds analogy principle. However, when the Prandtl number is distinctly different from unity, the spectral properties become different and new spectral laws are involved. In this respect, for the following developments, we shall rely upon the well known studies of [START_REF] Batchelor | Small scale variation of convected quantities like temperature in a 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] and Obukhov (1959), and reinterpreted afterwards by [START_REF] Warhaft | Passive scalars in turbulent flows[END_REF].

B.2

The case of small molecular Prandtl numbers

Each contribution k θsf s and k θ can be computed by integrating the spectrum E(κ) from their respective wave number ranges [κ c , κ e ] and [0, κ e ]. The first point is to calculate k θ using Eq. ( 41). As a result, it is a simple matter to obtain

k θ = ∞ 0 E θ (κ)dκ = C θ ǫ θ k C K ǫ ϑ H 0 E * (ϑ)dϑ = C θ ǫ θ k C K ǫ 1 -[1 + βϑ α H ] -γ (85)
so that the spectrum E θ (κ) given by Eq. ( 41) can be rewritten as

E θ (κ) = k θ k E(κ)H(κ H -κ) (1 -[1 + βϑ α H ] -γ ) (86)
In the first case where ϑ c < ϑ H , the ratio k θsf s /k θ can be computed easily using the spectrum given by Eq. ( 86) as k θsf s k θ =

B.3 The case of large molecular Prandtl numbers

The integration of E θ (κ) over the wave number range is given by

k θ = κ K 0 E θ 1 (κ)dκ + ∞ κ K E θ 2 (κ)dκ (90)
In the first case where ϑ c < ϑ K , then

k θsf s k θ = ϑ K ϑc E * ϑ 1 (ϑ)dϑ + ϑ S ϑ K E * ϑ 2 (ϑ)dϑ (
(ϑ K < ϑ C < ϑ S ) c ǫ θθ 1 (ϑ S < ϑ C ) (96)

C THE VISCOUS SUBLAYER LIMIT

The low Reynolds number extension of the model equations for the velocities and stresses have been defined in previous papers (Chaouat andSchiestel, 2005, 2013). In the present case, Eq. ( 54) for the variance dissipation-rate ǫ θ is modeled at low Reynolds number to approach walls using the same empirical approach that is of practical use for simulation of confined flows. As a result of modeling, the where Re sf s = k 2 sf s /(νǫ sf s ). This formulation at low Reynolds number as a function of the Prandtl number is inspired from previous papers, see for instance [START_REF] Shikazono | Second-moment closure for turbulent scalar transport at various Prandtl numbers[END_REF]; [START_REF] Nagano | Development of a two-equation heat transfer model based on direct simulations of turbulent flows with different Prandtl numbers[END_REF]; [START_REF] Kenjeres | Convective rolls and heat transfer in finite-length Rayleigh-Benard convection: A two-dimensional numerical study[END_REF]. Taylor series expansion in space for the fluctuating velocities u ′ i and fluctuating scalar variable θ ′ associated with the passive scalar show that the ratios ǫ/k sf s and ǫθ /k θsf s take finite values at the wall for x n = 0. In particular, the boundary condition for ǫ and ǫ θ at the wall are given by the relations ǫw = 2ν[∂ k

Figure 1 :

 1 Figure 1: Sketch of spectral splitting of kinetic energy of turbulence

As a result of

  modeling[START_REF] Schiestel | Towards a new partially integrated transport model for coarse grid and unsteady turbulent flow simulations[END_REF] Chaouat and Schiestel, 2005, 2012), the transport equation for the subfilter stress tensor (τ ij ) sf s can be written in the simple compact form as[START_REF] Chaouat | A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows[END_REF]

1 and c ǫ 2

 2 are the coefficients used in RANS with the values c ǫ 1 = 1.50 and c ǫ 2 = 1.90, respectively. This function introduced in Eq. (13) allows to sensitize the model to the filter width. Then, the coefficient c ǫ 2sf s can be considered as a dynamical parameter which draws the spectral distribution towards the prescribed equilibrium distribution. In other words, this term acts like a relaxation towards the Kolmogorov equilibrium spectrum[START_REF] Chaouat | Analytical insights into the partially integrated transport modeling method for hybrid Reynolds averaged Navier-Stokes equations-large eddy simulations of turbulent flows[END_REF]. The coefficient c ǫ 1sf s is the same as the one used in the corresponding RANS dissipation equation

Figure 2 :

 2 Figure 2: Sketch of spectral splitting of scalar variance

  as demonstrated in Eq. (81). Introducing the Kolmorogov scale η K = (ν 3 /ǫ)1/4 and the Prandtl number P r = ν/σ indicated in Refs.[START_REF] Batchelor | Small scale variation of convected quantities like temperature in a turbulent fluid. Part 1. General discussion and the case of small conductivity[END_REF][START_REF] Tennekes | A first course of turbulence[END_REF][START_REF] Hinze | Turbulence[END_REF][START_REF] Monin | Statistical Fluid Mechanics[END_REF], the scalar microscale can be computed by η θ = η K /P 3/4

Figure 3 :

 3 Figure 3: Spectral zones for small molecular Prandtl numbers where κ H = (ǫ/σ 3 ) 1/4 and L = k 3/2 /ǫ.

4 tFigure 4 : 2 r

 442 Figure 4: Spectral zones for large molecular Prandtl numbers where κ K = 1/η K and κ S = P 1/2 r /η K with η K = (ν 3 /ǫ) 1/4 and L = k 3/2 /ǫ.

a

  generalized gradient diffusion hypothesis (GGDH) first introduced by Daly and Harlow (1970) and often applied in practice by many authors such as[START_REF] Launder | On the computation of convective heat transfer in complex turbulent flows[END_REF],[START_REF] Suga | Improvement of second-moment closure for turbulent obstacle flow and heat transfer[END_REF],[START_REF] Kenjeres | LES,T-RANS and hybrid simulations of thermal convection at high Ra number[END_REF],[START_REF] Kenjeres | Capturing transient effects in turbulent flows over complex urban areas with passive pollutants[END_REF]

Figure 5 Figure 5 : 4 r 3 / 4 t 4 t 4 t

 5543444 Figure 5 displays the variation of the coefficient c ǫ θθ 2 sf s versus the location of the cutoff depending on

Figure 6 :

 6 Figure 6: Setup of the numerical channel flow simulations with various thermal boundary conditions. Both walls are heated.

  k 2 /(νǫ) is the turbulent Reynolds number, because the thermal scales are smaller than the dynamic Kolmogorov scales. At P r = 1, R 2 (P r ) = N DN S /N P IT M 2 = 451 and R 1 (P r ) = N DN S /N P IT M 1 = 113 while at P r = 10, R 2 (P r ) = N DN S /N P IT M 2 = 2378 and R 1 (P r ) = N DN S /N P IT M 1 = 594. This means

Figure 7 :Figure 8 :Figure 9 :Figure 10 :

 78910 Figure 7: Mean scalar field θ + = θ /θ τ in logarithmic coordinate versus the wall unit distance for various P r numbers. DNS : •; PITM1: . PITM2: . (a) P r = 0.1; (b) P r = 1; (c) P r = 10; R τ = 395.

Figure 11 :

 11 Figure 11: Contours of the instantaneous filtered variable Θ in the (x 1 , x 3 ) mid-plane illustrating the unsteady character of the scalar field. PITM2. Comparison can be made with DNS (Chaouat and Peyret, 2019).

  usual corresponding RANS equation but the coefficient c θθsf s 2 appearing in the destruction term is now a function of the ratio of the subfilter variance to the total variance k θsf s /k θ taking into account the cutoff wave number κ c and the turbulence length scale L = k 3/2 /ǫ. Considering different spectra E θ (κ) of the passive scalar depending on the Prandtl number range, the mathematical expression of the subfilter coefficient has been calculated as a function of the dimensionless variables involving the characteristic length-scales, the Prandtl number and the Peclet number. These foundation concepts in spectral space have then be used to develop the subfilter closure model extended to turbulent scalar fluctuations. In a first step, we have then considered for the turbulent field the second moment closure based on the transport equations of the subfilter scale stresses (τ ij ) sf s and the dissipation-rate ǫ, coupled with an algebraic diffusivity model for the scalar fluxes (τ iθ ) sf s used in the transport equations of the variance of the subfilter scalar field k θsf s and its dissipation-rate ǫ θ . Numerical simulations of fully turbulent channel flows have been performed on several grids for illustrating the capabilities of the present PITM model to reproduce passive scalar fields. Considering that a good prediction of turbulent scalar fluxes is a prerequisite for meaningful scalar variance study, thorough comparisons of heat fluxes with DNS data have been detailed with satisfactory results. The PITM yields these scalar field correlations in good agreement with DNS data for both grids. It has been verified that the subfilter contributions increase as the grid-size increases since larger scales must be modeled and vice versa, the resolved contributions decrease since less scales are computed. But as expected, the PITM model returns the total turbulent heat flux and scalar variance as free of any spectral cutoff location proving the entire self-consistency of the model. Moreover, some insights into the structural appearance of the scalar field interpreted as a marker for the fluid particles has been given. From a more general point of view, this pioneering work opens new routes of modeling for the simulation of turbulent flows including a passive scalar with a drastic reduction of the computational time and memory in term of number of grid points, in comparison with the demanding resources of highly resolved LES. Applications to practical and geophysical flows will be tackled in further works.A TRANSPORT EQUATION OF THE TWO POINT SCALARFLUCTUATION A.1 Transport equation of the two-point correlation tensor

Fig. 2 . 2 3

 22 Fig. 2. Using Eq. (37), it is then simple matter to compute the subfilter part of the variance k θsf s with the change of variable ϑ = κL by k θsf s k θ =

  case where ϑ c > ϑ H , then k θsf s /k θ = 0, but, in this case, there is no longer need to solve the ǫ θ equation ! The exact expression of the coefficient c ǫ θθ 2 sf s is thenc ǫ θθ 2 sf s =            c ǫ θθ 1 + ∆c ǫ θθ [1 + βϑ α c ] -γ -1 + βP e

c ǫ θθ 1 + ∆c ǫ θθ c θ ǫ θ k θ ν ǫ 1 /2 ln ϑ S ϑc 1 -

 111 the second case where ϑ K < ϑ c < ϑ S , then case where ϑ S < ϑ c , then k θsf s /k θ = 0. As a result, we obtain the exact expression for the coefficient c ǫ θθ 2 sf s valid in all casesc ǫ θθ 2 sf s =                          c ǫ θθ 1 + ∆c ǫ θθ [1 + βϑ α c ] -γ -1 + βRe

1 / 2

 12 sf s /∂x n ] 2 w and ǫθw = 2σ[∂ k θ 1/2 sf s /∂x n ] 2 w .The code developed by[START_REF] Chaouat | An efficient numerical method for RANS/LES turbulent simulations using subfilter scale stress transport equations[END_REF] is based on the finite volume technique so that the equations are solved in the grid cell around the grid-point of the mesh and not on the grid point itself.

Table 1 :

 1 Different levels of description in turbulence modeling including a passive scalar

		Level 1	Level 2	Level 3	Level 4
	Turbulence quantity	3D spectrum Spherical av-	Subfilter scale	Mean subfil-
			eraged spec-	part	ter scale part
			trum
	Variance of velocity
	fluctuations

  The subfilter tensorial dissipation rate ǫ ij is approached by 2/3ǫδ ij at high Reynolds number. Its transport equation reads[START_REF] Schiestel | Towards a new partially integrated transport model for coarse grid and unsteady turbulent flow simulations[END_REF] Chaouat and

	Schiestel, 2005, 2009, 2012)

Table 2 :

 2 Coefficients values used in the turbulence model.

and we choose c ǫ θθ 2 = 1.30. Note that there is no precise consensus in the scientific literature

[START_REF] Launder | Simulation of turbulence models and their applications[END_REF][START_REF] Yoshizawa | Statistical modelling of passive-scalar diffusion in turbulent shear flows[END_REF][START_REF] Nagano | A two-equation model for heat transport in wall turbulent shear flows[END_REF][START_REF] Ruffin | Investigation of characteristic scales in variable density turbulent jets using a second-order model[END_REF][START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF][START_REF] Hanjalic | Modelling turbulence in engineering and the environment. Second-moment route to closure[END_REF]

. With these numerical constants, the dimensionless number R takes on the value 0.75. The coefficients values used in the turbulence model are summarized in table 1. The

  so that, using the partial derivative operators defined in Appendix A[START_REF] Chaouat | From single-scale turbulence models to multiple-scale and subgridscale models by Fourier transform[END_REF]

		ǫ θ =	σ 4	∂ 2 φ θθ ∂X j ∂X j ξ=0 -σ	∂ 2 φ θθ ∂ξ j ∂ξ j ξ=0 = σ	∂θ ′ ∂x j	∂θ ′ ∂x j	(80)
	B INFLUENCE OF PRANDTL NUMBER IN THE SCALAR
		DISSIPATION EQUATION
	B.1	The case of molecular Prandtl numbers near unity
					2 j E ϑϑ	dκ		(79)
						ξ=0

  transport equation for ǫ θ reads∂ǫ θ ∂t = c ǫ θθ 1 sf s P θsf s ǫ θ k θsf s + c ǫ θk 1 sf s f ǫ θk 1 P sf s ǫ θ k sf s -c ǫ θk 2 sf s sf s /∂x n ] 2 , ǫθ = ǫ θ -2σ[∂ k θ 1/2sf s /∂x n ] 2 and f ǫ θk 1 is a damping function of the Prandtl and Reynolds numbers which has been calibrated as f ǫ θk 1 = 0.4(1-P r )/(1+P r ) exp[-( Re sf s /100)]

		ǫ θ k sf s ǫ	-c ǫ θθ 2 sf s	ǫ θ ǫθ k θsf s	+ J ǫ θ	(97)
	where ǫ = ǫ -2ν[∂ k	1/2