Contribution of direct numerical simulations to the budget and modelling of the transport equations for passive scalar turbulent fields with wall scalar fluctuations

Abstract

We perform direct numerical simulations (DNS) of the fully developed turbulent channel flow with a passive scalar subjected to constant time-averaged scalar fluxes at the wall. The first case (case I) considers a constant time-averaged scalar flux $\left \langle q_w \right \rangle$ along with the specific condition of a non-fluctuating scalar $\theta _w$ imposed at the wall implying zero wall scalar fluctuation $\theta '_w=0$ and zero variance (isoscalar boundary condition) whereas the second case (case II) accounts for a constant instantaneous scalar flux $q_w$ in time and space leading to zero gradient of scalar fluctuation along the normal to the wall $(\partial \theta '/ \partial x_n)_w = 0$ (isoflux boundary condition) implying a non-zero variance. The friction Reynolds number takes on the value $R_\tau =395$ and the molecular Prandtl numbers considered are $P_r=0.01$, 0.1, 1 and 10. The purpose is to investigate the effect of the wall scalar fluctuations on the scalar field. Emphasis is put on the mean passive scalar $\left \langle \theta \right \rangle$, the half-scalar variance $k_\theta$, the turbulent scalar fluxes $\tau _{i \theta }$, the correlation coefficients $R_{i \theta }$, the passive scalar to dynamic time-scale ratio $\mathcal {R}$, the turbulent Prandtl number $Pr_t$ and higher-order scalar statistics. Systematic comparisons between these two scalar fields are undertaken. As a result of interest, it is found that the mean scalar remains almost the same whatever the type of the boundary layer condition but not the scalar variance. The budgets of transport equations for the half-scalar variance and turbulent fluxes reveal that some scalar quantities such as the dissipation-rate and the molecular diffusion terms are highly modified in the near-wall region but not really the production, diffusion and the scalar-pressure gradient correlation terms in a first approximation. Visualization of the instantaneous scalar fields indicates that the topology of the structures is strongly modified as the Prandtl number increases from $P_r=0.01$ up to $10$. Finally, it is shown how to use this DNS database to devise and calibrate scalar flux equation models in the framework of both first and second moment closures. This study suggests that accounting for wall scalar fluctuations should be considered in the simulation of turbulent flows involving fluid and solid combinations at the interface and provides a useful high resolution DNS database.

1. Introduction

Turbulent flows accounting for the transport of a passive scalar are often encountered in nature involving for instance the pollution dispersal in the atmosphere that has been of concern in the environment and climate change but also in industrial processes including alumina refinery plant, combined cycle power plant, heat exchangers, mixing and other systems. The passive scalar hypothesis holds when the turbulent velocity field mainly governs the transport of the scalar field while the mutual influence of the scalar field on the velocity field is rather weak and can be neglected in a first approximation. In various flow configurations involving heat transfer between impermeable walls, it is of importance to accurately determine the thermal boundary condition on the wall resulting from the interaction between the heat conduction in the solid and the turbulent flow motion that appreciably affects the flow solution. This problem is particularly acute in industrial applications such as for instance liquid metal-cooled reactors when the temperature fluctuations at the wall may lead to thermal fatigue failure of solid structures. This problem may appear in the mixing of two flows at different temperatures merging at a T-junction (Georgiou & Papalexandris 2018). For instance, among others, thermal striping can arise in certain liquid metal-cooled reactors. As the temperature fluctuations play an important role in industrial devices, these effects must be deeply investigated and calibrated. As confirmed by the experimental study of the wall temperature fluctuations conducted by Mosyak, Pogrebnyak & Hetsroni (2001), the flow solution at the interface for the solid–fluid conjugate system depends on the thermal effusivity ratio of the gas to the solid given by $K=a/a_w$, where the thermal effusivity is defined as $a=\sqrt {\rho c_p \kappa }$ involving the density $\rho$, the specific heat at constant pressure $c_p$ and the scalar conductivity $\kappa$, the subscript $w$ referring to the properties of the wall. Usually, the ratio of the effusivity of the gas to that of structural materials is very small $K \leq 10^{-3}$ but not for liquids, where it is of order unity (Kasagi, Kuroda & Hirata 1989). For a turbulent flow subjected to heat fluxes through the wall, two different boundary conditions can be usually considered depending on the fluid and solid properties. The first one is the isoscalar boundary condition with a constant wall temperature (case I) implying that its fluctuation at the wall $\theta '_w$ is zero and corresponds to the usual case where $K \ll 1$, i.e. $a \ll a_w$ with a large wall thickness at the interface. The second is the isoflux boundary condition implying non-zero temperature fluctuations at the wall (case II) and corresponds to the case where $K \gg 1$, i.e. $a \gg a_w$ with an infinitesimal wall thickness. In this framework, Kasagi et al. (1989) have developed an unsteady streamwise pseudo-vortical motion model to investigate the near-wall behaviour of the statistical quantities of temperature fluctuations for various Prandtl numbers with the influence of the thermal properties and thickness of the wall. The numerical results returned by this model were found to be in agreement with the experimental study of the wall temperature fluctuations (Mosyak et al. 2001), indicating that the fluctuating temperatures near the wall are of higher intensity for the isoflux condition (case II) than for the isothermal condition (case I).

Direct numerical simulation (DNS) development began with the pioneering work of the Stanford team with the investigation of the fully developed turbulent flow in a plane channel (Kim, Moin & Moser 1987) for the Reynolds number $R_\tau =u_\tau \delta /\nu =180$ based on the friction velocity $u_\tau$, half-channel width $\delta$ and dynamic viscosity $\nu$, and is now a very valuable tool (Lee & Moser 2015; Paul, Papadakis & Vassilicos 2018) thanks to the rapid increase of super-computer power as well as the development of computational techniques including vectorization and parallelization (Schiestel & Chaouat 2022). Numerical simulations of turbulent flows are useful also in this framework of heat transfer for determining the thermal flow characteristics. It is the best tool to apply with the highest accuracy possible, considering that measurements of concentration of passive contaminants are very difficult to perform. The case of the channel flow heated on both walls with a constant imposed heat flux $q_w$ or subjected to a wall temperature difference $\Delta T_w$ and zero wall temperature fluctuation was first performed by Kim & Moin (1989) and was since deeply investigated by several authors where a non-exhaustive list is summarized in table 1. These papers provided useful information about the mean temperature, the variance, the turbulent heat fluxes and other quantities. Among these, the two recent papers of Alcántara-Ávila, Hoyas & Pérez-Quiles (2018) and Alcántara-Ávila & Hoyas (2021) concern simulations for the Reynolds numbers $R_\tau =500$, 1000 and 2000 and several Prandtl numbers varying from low, medium to high values. The budget of the temperature variance and its dissipation rate were obtained in addition to the effect of the Reynolds number on the statistics. In particular, it was shown that the maximum of the intensity of the thermal field does not increase with the Reynolds number for the highest Prandtl numbers. This well-documented data base considers only one type of boundary condition, the present work extends the analysis to two differing boundary conditions including both with and without wall temperature fluctuations (or wall scalar fluctuations). Some authors reported studies with wall temperature fluctuations but these papers do not contain comprehensive data (Lu & Hetsroni 1995; Tiselj et al. 2001; Flageul et al. 2015). These authors showed that the effect of the wall scalar fluctuations is limited to the wall region. For these simulations, the budget of temperature variance and heat fluxes were performed for the Reynolds and Prandtl numbers $R_\tau =171$, $P_r=1$; $R_\tau =180$, $P_r=0.72$; $R_\tau =150$, $P_r=0.71$, respectively, and have indicated that the wall temperature fluctuations increase the scalar variance at the wall and reduce the magnitude of the molecular diffusion and dissipation of the scalar variance. This case was investigated more recently for the Reynolds number $R_\tau =395$ and the Prandtl numbers $P_r=0.01$, 0.1, 1, 10 (Chaouat 2018; Chaouat & Peyret 2019). The results have revealed how the thermal characteristics of the flow evolve as the molecular Prandtl number varies from low to high values but the budgets of the scalar transport equations, the passive scalar to dynamic time-scale ratio $\mathcal {R}$, amongst others, were not performed.

Table 1. The DNS and under-resolved DNS of turbulent channel flow subjected to different boundary conditions, constant heat flux $q_w$ or wall temperature difference $\Delta T_w$ with zero wall temperature fluctuation $T'_w=0$ performed for various Reynolds and Prandtl numbers.

The objective of this paper is twofold. On the one hand, it is to thoroughly investigate the effect of the wall scalar fluctuations on the scalar field with emphasis on the budget data of the transport equations for the half-scalar variance $k_\theta$ and the turbulent scalar fluxes $\tau _{i \theta }$ for the Reynolds number $R_\tau =395$ and the Prandtl numbers $P_r=0.01$, 0.1, 1 and 10, and on the other hand, to provide a useful high resolution DNS database for users and researchers aiming to validate turbulence models with transfer of a passive scalar. This work will contribute for a better understanding of the scalar fluctuation effect for $R_\tau =395$ and very low to high Prandtl numbers ranging from 0.01 to 10 in a descriptive framework which has not yet been reported in the open literature. Physically, the range of Prandtl numbers has been chosen to be representative of heat transfers occurring in industrial applications involving liquid or gases. In this framework, DNSs of the turbulent channel flow subjected to constant scalar fluxes as shown in figure 1 are performed on several grids of very high resolution. In addition, it will be shown how to use this DNS database to devise and calibrate scalar flux equation models. This database will constitute a useful resource for researchers involved in Reynolds averaged Navier–Stokes model equations (RANS) (Hamba 2004; Schiestel 2008; Gatski 2009; Hanjalic & Launder 2011), large eddy simulation (LES) (Lesieur & Métais 1996; Meneveau & Katz 2000; Abe & Suga 2001) and also hybrid RANS/LES methodologies (Chaouat 2017b; Chaouat & Schiestel 2005, 2012, 2021a,b), considering that these methodologies are complementary tools from each other, as mentioned recently by Schiestel & Chaouat (2022).

Figure 1. Set-up of the numerical channel flow simulations subjected to constant scalar fluxes from both walls in the normal direction $x_n$. Two different scalar boundary conditions are applied leading to two separate cases for the fluctuating scalar variance (a) $\theta '_{w}=0$ (case I); (b) $q'_w= \kappa (\partial \theta '/ \partial x_n)_w =0$ (case II).

2. Governing equations and boundary conditions

2.1. Equations

The governing equations for the velocity and scalar fields are the incompressible three-dimensional mass conservation, Navier–Stokes equations and the transport equation for the passive scalar. As usual, the fully developed turbulent channel flow is driven by a constant streamwise pressure gradient term and the passive scalar is subjected to a constant scalar flux at the lower and upper walls. The Navier–Stokes equations for the velocity field expressing the conservation of mass and momentum are

(2.1)\begin{equation} \frac{\partial u^+_j } {\partial x^*_j } = 0 \end{equation}

and

(2.2)\begin{equation} \frac{\partial u^+_i } {\partial t^*} + \frac{\partial} {\partial x^*_j} (u^+_i u^+_j) ={-}\frac{\partial p^+} {\partial x^*_i} + \frac{1} {R_\tau } \frac{\partial ^2 u^+_i } {\partial x^*_j \partial x^*_j} +G_i, \end{equation}

where, in these equations, $u_i$ denotes the velocity and $p$ is the pressure. The time $t$, coordinate $x_i$ and flow variables are normalized using the channel half-width $\delta$, the friction velocity $u_\tau$ and the kinematic viscosity $\nu$ such as $t^* = t u_\tau /\delta$, $x^*_i = x_i/\delta$, $u^+_i= u_i/u_\tau$, $p^+=p/\rho u^2_\tau$. The quantity $G_1$ included in this latter equation denotes the mean pressure gradient term necessary to balance friction at the upper and lower walls allowing us to get a periodic condition between the inlet and outlet sections of the channel and takes the value unity in dimensionless form. Because of the scalar flux condition, the mean scalar variable $\langle \varTheta \rangle$ as well as the bulk mean scalar variable $\langle \varTheta _m \rangle$, where the angular bracket $\langle.\rangle$ denotes the averaging in time and space in the homogeneous directions, increase linearly in the $x_1$ direction. With the aim to keep a constant mean value in the channel along the streamwise direction, as proposed by Kawamura et al. (1998, 1999), a change of variable is made by introducing the new scalar variable $\theta$ such that

(2.3)\begin{equation} \theta = \frac{\partial \left \langle \varTheta_m \right \rangle } {\partial x_1} x_1 - \varTheta, \end{equation}

where the gradient ${\partial \langle \varTheta _m \rangle }/ {\partial x_1}$ takes on a constant value. With the above transformation (2.3), the transport equation for the dimensionless passive scalar $\theta ^+$ reads

(2.4)\begin{equation} \frac{\partial \theta^+ } {\partial t^*} + \frac{\partial} {\partial x^*_j} (\theta^+ u^+_j) = \frac{1} {R_\tau P_r} \frac{\partial ^2 \theta^+ } {\partial x^*_j \partial x^*_j} +Q, \end{equation}

where, in this equation, the variable $\theta$ is normalized by the wall characteristic transfer of the passive scalar (or friction temperature) such that $\theta ^+ = \theta /\theta _\tau$ with $\theta _\tau = q_w/(\rho c_p u_\tau )$, the velocity $u^+_i$ is given by (2.2) and $Q$ denotes the source term given in dimensionless form by $Q= {u^+_1} /{U^+_b} = {u^+_1} \partial \langle \varTheta ^+_m \rangle / {\partial x^*_1}$, $U^+_b=U_b/u_\tau$ where $U_b=(1/2\delta ) \int _{0}^{2 \delta } u_1 \,{\rm d}x_3$ is the average velocity over the channel cross-section. Its physical meaning corresponds to the mean scalar gradient necessary to balance wall scalar fluxes. The scalar flux applied at the walls is given by $q_w=-\kappa (\partial \theta /\partial x_3)_w$, where $\kappa = \rho c_p \nu /P_r$ stands for the scalar conductivity while the scalar diffusivity is given by $\sigma =\kappa /(\rho c_p)=\nu /P_r$. The averaged transport equation for the passive scalar in a steady flow regime reads

(2.5)\begin{equation} 0 = \frac{1} {R_\tau P_r} \frac{\partial ^2 \left\langle\theta^+\right\rangle } {\partial x^*_j \partial x^*_j} -\frac{\partial \tau^+_{j \theta} } {\partial x^*_j} + \left \langle Q \right \rangle, \end{equation}

where, in this equation, the turbulent scalar fluxes $\tau ^+_{i \theta } = \langle u'^{+}_i \theta '^+ \rangle$ play a crucial role in flow behaviour involving transfer of the passive scalar. The equation governing the wall-normal turbulent scalar flux is obtained by integrating (2.5) over the wall distance from $0$ to $x^+_3 = x_3 u_\tau /\nu$ leading to

(2.6)\begin{equation} \tau^+_{3 \theta} ({x^+_3}) = \frac{1} { P_r} \frac{\partial \left\langle\theta^+\right\rangle} {\partial x^+_3 } - 1 + \frac{1} {R_\tau U^+_b} \int_{0}^{x^+_3} \left \langle u^+_1 \right \rangle {\rm d}x'^+_3 . \end{equation}

2.2. Boundary conditions

The velocity boundary conditions at the lower and upper walls for $x_3=0$ and $2\delta$ are no-slip velocity conditions $u^+_i=0$. Constant scalar fluxes $q_w$ are applied on the lower and upper walls as shown in figure 1. In the general case of a solid–fluid interface at $x_3=0$, ($x_3 < 0$ for the solid, and $x_3 > 0$ for the fluid), it can be recalled that the dimensionless equation to be solved in the solid reads

(2.7)\begin{equation} \frac{\partial \theta^+ } {\partial t^*} = \frac{1} {R_\tau P_r} \frac{\partial ^2 \theta^+ } {\partial x^{**}_j \partial x^{**}_j} , \end{equation}

where $x^{**}_i = (\sigma /\sigma _w ) x^*_i$. The matching conditions for the scalar and the normal scalar flux at the interface are

(2.8)\begin{equation} (\theta^+ )_{x^*_3=0^+} = (\theta^+ )_{x^*_3=0^-} \end{equation}

and

(2.9)\begin{equation} K \left(\frac{\partial \theta^+} {\partial x^*_3 } \right)_{x^*_3=0^+}= \left(\frac{\partial \theta^+} {\partial x^{**}_3 } \right)_{x^*_3=0^-} \end{equation}

(Luikov 1968). As discussed in detail by Kasagi et al. (1989), the two major parameters for the conjugate heat transfer system are the thermal effusivity ratio $K=\sqrt {\rho c_p \kappa /\rho _w c_{{p}_w} \kappa _w}$ and the dimensionless wall thickness $d^{**}$ of the solid. In the present case, a simple way consists in considering the two limiting values of the thermal effusivity ratio (Kasagi et al. 1989; Li et al. 2009; Chaouat & Peyret 2019). Firstly, when $K$ goes to zero, i.e. $\theta _w$ is close to a constant value so that $\theta '_w=0$ at the wall leading therefore to zero root-mean-square (r.m.s.) fluctuations $\left \langle \theta ' \theta '\right \rangle _w=0$ (case I). Secondly, when $K$ goes to infinity, i.e. $q_w$ is imposed to a constant value in time and space implying that $q'_w=-\kappa (\partial \theta /\partial x_3)_w=0$, leading to non-zero r.m.s. fluctuations $\left \langle \theta ' \theta ' \right \rangle _{w} \ne 0$ (case II). Consequently, the r.m.s. intensity of the scalar fluctuations also verifies the condition ($\partial \langle \theta ' \theta ' \rangle /\partial x_3)_{w}=0$. This condition means that the variance of the passive scalar is almost constant along the normal to the wall. Mathematically, the isoscalar and isoflux boundary conditions correspond to the Dirichlet and Neumann boundary conditions, respectively, and sometimes require numerical procedures to be applied at the solid–fluid interface (Errera & Chemin 2013; Flageul et al. 2015).

2.3. Numerical procedure

For all simulations, the dimension of the channel in the streamwise, spanwise and normal directions along the Cartesian axes $x_1$, $x_2$, $x_3$ are $L_1= 6.4 \delta$, $L_2=3.2 \delta$ and $L_3=2 \delta$. It has been checked that the computational box is sufficiently long to allow the two-point correlation tensor of the fluctuating velocities in the streamwise direction to return to zero from the inlet to the outlet sections of the channel. The number of grid points $N_i$, the spacings $\varDelta _i$ as well as the Batchelor length scale $\eta _\theta$ relative to the Kolmogorov length scale $\eta _\kappa = (\nu ^3/\epsilon )^{1/4}$ are given in table 2. The grid resolution is determined in order to solve both the Kolmogorov scale $\eta _\kappa$ and the Batchelor length scale $\eta _\theta$ (Batchelor 1971; Tennekes & Lumley 1972) which approaches $\eta _\kappa$ when $P_r$ is of the order of unity, $\eta _\theta =(\sigma ^3/\epsilon )^{1/4}=\eta _\kappa /P^{3/4}_r$ at small Prandtl numbers and $\eta _\theta = (\nu \sigma ^2 /\epsilon )^{1/4}=\eta _\kappa /P^{1/2}_r$ at large Prandtl numbers. In particular, at $P_r=0.01$, $\eta _\theta \approx 31.6 \, \eta _\kappa$, at $P_r=0.1$, $\eta _\theta \approx 5.62 \, \eta _\kappa$, and at $P_r=1$, $\eta _\theta \approx \eta _\kappa$ but at $P_r=10$, $\eta _\theta \approx 0.316 \eta _\kappa$. The grid resolution is low at small Prandtl numbers but high at large Prandtl numbers according to the power law of the Prandtl number. It is then a simple matter to see that the number of grid points varies according to the law $N \propto Re^{9/4}_t P^{3/2}_r$ at large Prandtl numbers where $Re_{t}= k^2/(\nu \epsilon )$ is the turbulent Reynolds number. The computational time can be estimated if we assume that the turbulent Reynolds number is proportional to the bulk Reynolds number leading to $t \propto Re^{11/4}_{t} P^{3/2}_r$. Thus, the DNS of scalar turbulence becomes much more stringent to perform as the Prandtl number increases. The literature indicates that most of direct numerical simulations have been limited to Prandtl number smaller than 10 for $R_\tau \geqslant 395$. For the Reynolds number and Prandtl number values studied here, the grid numbers of the meshes are $M_1$ with the resolution $512 \times 256 \times 256$ for $P_r=0.01$, $0.1$ and 1, and $M_2$ with the resolution $1024 \times 512 \times 1024$ for $P_r=10$. The grid spacings are given by $\varDelta ^+_i = R_\tau L_i /N_i \delta$ leading to $\varDelta ^+_1=\varDelta ^+_2 \approx 5$ for $M_1$ and $\varDelta ^+_1=\varDelta ^+_2 \approx 2.5$ for $M_2$. For these meshes, $\varDelta ^+_3$ in the normal direction to the wall $x_3$ verifies the constraint condition $\varDelta ^+_3 < 2$ for $M_1$ and $\varDelta ^+_3 < 0.1$ for $M_2$, this value being very small to accurately compute the thermal boundary layer.

Table 2. Domain size, grid points, grid resolution and Batchelor length scale $\eta _\theta$ relative to the Kolmogorov length scale $\eta _\kappa$ for the Prandtl numbers $P_r=0.01$, $P_r=0.1$, $1$ and $10$; $N_1$, $N_2$ and $N_3$ are the number of grid points in the streamwise, spanwise and wall-normal directions, respectively, for the meshes $M_i$; DNS is performed at the friction Reynolds number $R_\tau = u_\tau \delta /\nu =395$ corresponding to the bulk Reynolds number $Re_m= 2 U_b \delta /\nu =13\,750$.

A grid sensitivity study has been conducted to ensure that the grid resolution is adequate. In particular, in regard to DNS performed at the Reynolds number $R_\tau =395$ for the Prandtl number $P_r=10$ (Chaouat 2018; Chaouat & Peyret 2019), the grid in the present case is more refined to compute the budgets of the transport equations for $k_\theta$ and $\tau _{i \theta }$ involving higher-order statistics. Hence, the additional Mesh $M_1$ with an extremely high resolution $1024 \times 512 \times 1024$, especially in the normal direction to the wall, has been used to capture adequately the correlation of the fluctuating scalar gradients $\langle ( \partial \theta ' / \partial x_i ) ( \partial \theta ' / \partial x_i ) \rangle$ or even the correlation of the fluctuating velocity-scalar gradients $\langle ( \partial u'_i / \partial x_i ) ( \partial \theta ' / \partial x_i ) \rangle$ associated with the fine structures that evolve locally in time and space. Note that the ratio of the number of grid points in the normal direction for the meshes $M_1$ and $M_2$ is set to $(N_3)_{M_1}/(N_3)_{M_2}=4$ to take into account the decrease in the Batchelor scale as the Prandtl number increases. Finally, it has been checked also that higher resolution of the grid size in all directions for each mesh yields a change in the fluctuation intensity less than the order of a few per cent.

The governing equations (2.2) and (2.4) are integrated using the explicit Runge–Kutta scheme with fourth-order accuracy in time and are solved in space by means of a conservative centred fourth-order scheme in the $x_1$, $x_2$, $x_3$ directions as described in Appendix A. The time step is adjusted to get a constant Courant–Friedrichs–Lewy number of unity. This combined time–space numerical scheme of high-order accuracy has shown very good numerical properties in performing DNS of turbulence with a passive scalar (Chaouat 2018; Chaouat & Peyret 2019). The computational fluid dynamics (CFD) code developed by Chaouat (2011) is based on the finite volume technique where the mean variables are evaluated at the centre of the computational cell whereas the convective and diffusive fluxes are computed at the interfaces surrounding the cell. The numerical code is highly optimized with message passing interface (MPI) thanks to the compact and explicit formulation of all involved stencils. The present DNSs were run for a sufficiently long time to get fully converged statistics.

3. Transport equations for scalar variables

The transport equations for the scalar quantities such as the half-variance $k_\theta$ and the turbulent fluxes $\tau _{i \theta }$ are written in the fully developed turbulent flow and scalar field for a steady statistical flow regime in a non-dimensionalized form using the wall unit $x^+_i$. The estimate of the terms appearing in equations (3.1) and (3.2) is given for quasi-isotropic turbulence in Appendix B for the sake of clarity.

3.1. Transport equation for the half-scalar variance

The transport equation for the half-scalar variance $k^+_\theta = \langle \theta '^+\theta '^+ \rangle /2$ reads

(3.1) \begin{align} 0 &= \underbrace{ \tau^+_{1 \theta} \frac{\partial \langle\varTheta^+_m \rangle } {\partial x^+_1 } - \tau^+_{j \theta} \frac{\partial \langle\theta^+ \rangle } {\partial x^+_j } }_{P^+_\theta= O(1)} - \underbrace{ \frac{1} {2} \frac{\partial} {\partial x^+_j} \langle u'^+_j \theta'^{{+}2} \rangle }_{T^+_\theta= O(1)} + \underbrace{ \frac{1} {P_r} \frac{ \partial^2 k^+_\theta} {\partial x^{2+}_j} }_{d^+_{\theta}= O ({Re}^{{-}1}_t P^{{-}1}_r )} \nonumber\\ &\quad - \underbrace{ \frac{1} {P_r} \left \langle \frac{\partial \theta'^+} {\partial x^+_j} \frac{\partial \theta'^+} {\partial x^+_j} \right \rangle}_{\epsilon^+_\theta = O(1)}, \end{align}

where, in this equation, $P^+_\theta$ denotes the production, in which the first term results from the linear increase of the bulk scalar variable, while the second term is caused by the interaction between the turbulent heat flux and the mean scalar gradient, $T^+_\theta$ is the turbulent diffusion due to the correlation of the scalar velocity, $d^+_\theta$ is the molecular diffusion and $\epsilon ^+_\theta$ is the dissipation rate of the half-scalar variance. Note that the first production term is negligibly small compared with the second term because ${\partial \left \langle \varTheta ^+_m \right \rangle }/{\partial x^+_1 } \ll {\partial \left \langle \theta ^+ \right \rangle } /{\partial x^+_3 }$. As expected, the molecular diffusion term $d^+_{\theta }$ varies according to the Reynolds number $Re^{-1}_t$ so that it is negligible at low Reynolds numbers in comparison with the other terms, in particular, far away from the wall. The transport equation for the dissipation rate $\epsilon ^+_\theta$ appearing in the last term in (3.1) is given in Appendix C but not discussed in detail here for the sake of conciseness. In addition, the transport equation for the turbulent kinetic energy $k$ is recalled in Appendix D for comparison purpose with $k_\theta$.

3.2. Transport equation for the turbulent scalar flux

The transport equation for the turbulent scalar flux $\tau ^+_{i \theta } = \langle u'^+_i \theta '^+ \rangle$ reads

(3.2)\begin{align} 0 &= \underbrace{ \tau^+_{1 i} \frac{\partial \left\langle\varTheta^+_m \right\rangle } {\partial x^+_1 } }_{P^{+(1)}_{i \theta}=O(1)} \underbrace{ - \tau^+_{i j } \frac{ \partial \left \langle \theta^+ \right \rangle } {\partial x^+_j } }_{P^{+(2)}_{i \theta}=O(1)} \underbrace{ - \tau^+_{j \theta} \frac{ \partial \left \langle u^+_i \right \rangle } {\partial x^+_j } }_{P^{+(3)}_{i \theta}=O(1)} \underbrace{ - \frac{\partial} {\partial x^+_j} \left \langle u'^+_iu'^+_j \theta'^+ \right \rangle }_{T^+_{i \theta}= O(1)} \nonumber\\ &\quad + \underbrace{ \frac{\partial} {\partial x_j} \left \langle \theta'^+ \frac{\partial u'^+_i} {\partial x^+_j} \right \rangle} _{d^{+(1)}_{i \theta}=O \left ({Re}^{{-}1/2}_t \right )} + \underbrace{ \frac{1} {P_r} \frac{\partial} {\partial x_j} \left \langle u'^+_i \frac{\partial \theta'^+} {\partial x^+_j} \right \rangle} _{d^{+(2)}_{i \theta}= O \left ( {Re}^{{-}1/2}_t P^{{-}1/2}_r \right) } \underbrace{ - \left \langle \theta'^+ \frac{ \partial p'^+} {\partial x^+_i} \right \rangle }_{\varPi^+_{i \theta}=O(1)} \nonumber\\ &\quad - \underbrace{ \left (1+\frac{1} {P_r} \right) \left \langle \frac{\partial u'^+_i } {\partial x^+_j } \frac{\partial \theta'^+ } {\partial x^+_j } \right \rangle }_{\epsilon^+_{i \theta} = O \left ( (P^{{-}1/2}_r +P^{1/2}_r ){Re}^{{-}1/2}_t \right )}, \end{align}

where, in this equation, $\tau ^+_{ij} = \left \langle u'^+_i u'^+_j \right \rangle$, $P^{+(1)}_{i \theta }$, $P^{+(2)}_{i \theta }$ and $P^{+(3)}_{i \theta }$ denote the production terms caused by the streamwise bulk scalar gradient, the interaction between the turbulent scalar flux with the mean velocity or scalar gradients, $T^+_{i \theta }$ is the turbulent diffusion caused by the correlation of the velocity-scalar variable, $d^{+(1)}_{i \theta }$ and $d^{+(2)}_{i \theta }$ are the viscous and scalar molecular diffusions, $\varPi ^+_{i \theta }$ is the correlation of the scalar-pressure gradients, $\epsilon ^+_{i \theta }$ denotes the scalar dissipation term. Note that the scalar-pressure gradient correlation term $\varPi ^+_{i \theta }$ is usually split into a pressure scalar gradient correlation $\varPhi ^+_{i \theta }$ and a pressure diffusion term $d^{+(3)}_{i \theta }$ as

(3.3)\begin{equation} \underbrace{ - \left \langle \theta'^+ \frac{ \partial p'^+} {\partial x^+_i} \right \rangle }_{\varPi^+_{i \theta} } = \underbrace{ \left \langle p'^+ \frac{ \partial \theta'^+} {\partial x^+_i} \right \rangle }_{\varPhi^+_{i \theta} } \underbrace{ - \left \langle \frac{\partial ( \theta'^+ p'^+) } {\partial x^+_i} \right \rangle }_{d^{+(3)}_{i \theta}} , \end{equation}

where the first term $\varPhi ^+_{i \theta }$ is the well-known pressure-scalar gradient correlation term leading to a better insight into the acting mechanisms involving the pressure fluctuation whereas the second term $d^{+(3)}_{i \theta }$ is usually interpreted as the turbulent diffusion due to the pressure fluctuations. Once again, the diffusion terms $d^{+(1)}_{i \theta }$ and $d^{+(2)}_{i \theta }$ as well as the scalar dissipation terms varying as ${Re}^{-1/2}_t$ are negligible at high turbulent Reynolds numbers. The order of magnitude of $\epsilon ^+_{i \theta }$ has been estimated considering that this term reduces to zero for isotropic turbulence in contrast with the dissipation of the scalar variance.A wall asymptotic analysis of the several terms appearing in these equations is made in Appendix E.

4. Numerical results and discussion

4.1. Kolmogorov and Batchelor turbulence scales

Figure 2 displays the evolution of the Kolmogorov and Batchelor dimensionless scales $\eta ^+_\kappa$ and $\eta ^+_\theta$, respectively, for Prandtl numbers ranging from 0.01 to 10 in addition to the grid spacing $\varDelta ^+_3$ for the meshes $M_1$ and $M_2$. It is found that these scales $\eta _\kappa$ and $\eta _\theta$ are approximately constant and minimum in the near-wall region but increase in the outer region when moving away from the wall due to the fact that the turbulence field becomes more isotropic. As expected, the grid spacing $\varDelta ^+_3(M_2)$ given in table 2 is smaller than $\eta ^+_\theta$ for $P_r=10$, $\varDelta ^+(M_1) < \eta ^+_\theta$ in the near-wall region but is, however, of the same order of magnitude as $\eta ^+_\theta$ in the central region for $P_r=1$, and $\varDelta ^+(M_1) < \eta ^+_\theta$ for $P_r=0.1$ and 0.01, confirming that the grid resolution is appropriate for performing the DNSs. Formally, the Batchelor length scale $\eta _\theta$ considered as a function of the Kolmogorov length scale $\eta _\kappa$ and the molecular Prandtl number $P_r$ (see § 2.3) is not modified by the wall boundary condition because the dynamic of the fluid remains the same for a passive scalar.

Figure 2. Kolmogorov length scale $\eta ^+_\kappa$ and Batchelor length scale $\eta ^+_\theta$ computed for the Prandtl numbers $P_r= 0.01$, 0.1, 1 and 10 in logarithmic coordinates vs the wall unit distance: - - -, $\varDelta ^+_3(M_1)$; —, $\varDelta ^+_3(M_2)$; $\blacktriangledown$ (cyan), $\eta ^+_\theta$, $(P_r=0.01)$; $\blacktriangle$ (blue), $\eta ^+_\theta$, $(P_r=0.1)$; $\blacksquare$ (green), $\eta ^+_\theta = \eta ^+_\kappa$, $(P_r=1)$; $\bullet$ (red), $\eta ^+_\kappa$, $(P_r=10)$. Here, $R_\tau =395$.

4.2. Statistics

The statistics for the mean scalar $\theta$, the half-scalar variance, $k_\theta$, the scalar fluxes $\tau _{i \theta }$, the correlation tensor $R_{i \theta }$ and the budget equations for $k_\theta$, $\tau _{1 \theta }$ and $\tau _{3 \theta }$ are investigated for the Prandtl numbers $P_r =0.01$, 0.1, 1 and 10. As mentioned above, DNSs are performed considering both the isoscalar and isoflux wall boundary conditions to enable comparisons between these two respective scalar fields. For sake of clarity, $\theta$ can be interpreted as the temperature, and its fluctuation is denoted $\theta '$ but it designates here any passive scalar variable associated with concentration or mass transfer transported by the turbulent flow. In the following, the turbulence and scalar quantities are normalized using the wall unit $x^+_i=x_i u_\tau /\nu$.

4.3. Mean scalar field

Figure 3 displays the mean scalar profile $\left \langle \theta \right \rangle$ in logarithmic coordinates vs the wall unit distance for the Prandtl numbers $P_r=0.01$, 0.1, 1 and 10. The conduction region penetrates less deeply into the core region of the channel when the Prandtl number increases because of the decrease of the scalar diffusivity $\sigma = \nu /P_r$, giving rise to an increase of $\left \langle \theta \right \rangle$. The logarithmic region begins to appear at $P_r=1$ and is highly extended for $P_r=10$. It is observed that the mean scalar profile remains the same whatever the type of the boundary condition applied at the wall, isoscalar or isoflux boundary conditions, indicating that there is no effect of the fluctuations at the wall, as also suggested by Kasagi et al. (1989) for an open channel flow and Chaouat & Peyret (2019) for the bounded channel flow, according to (E10). Figure 3(c) shows in addition the mean velocity $\langle u_1 \rangle$. As expected, the two curves plotted for $\langle \theta \rangle$ for $P_r=1$ and the one for $\langle u_1 \rangle$ are exactly the same, emphasizing the analogy between these two respective fields. Obviously, the mean velocity profile remains exactly the same at all Prandtl numbers because the scalar is passive but it considerably differs from the scalar profiles themselves for the other Prandtl numbers.

Figure 3. Mean scalar field $\left \langle \theta \right \rangle$ and mean velocity $\left \langle u_1 \right \rangle$ plotted only in case (c) for $P_r=1$ in logarithmic coordinate vs the wall unit distance. Here, $\theta '_w=0$ (case I), $\bullet$ (maroon); $q'_w=0$ (case II), $\blacklozenge$ (blue); $\langle u_1 \rangle$, $\blacktriangle$ (cyan). Panels show (a) $P_r=0.01$; (b) $P_r=0.1$; (c) $P_r=1$; (d) $P_r=10$; $R_\tau =395$.

4.4. The r.m.s. of the scalar variance

Figure 4 exhibits the r.m.s. of the scalar variance for the Prandtl numbers $P_r=0.01$, 0.1, 1 and 10, showing clearly the effect of the wall fluctuations on the scalar field in the immediate vicinity of the wall. A different scale in the axis depending on the Prandtl number has been chosen here in order to highlight the differences of behaviour near the wall, the smaller the Prandtl number is, the broader is the near-wall peak. As expected, the r.m.s. intensity reduces to zero at the wall for the isoscalar boundary condition (case I) but not at all for the isoflux boundary condition (case II) where it is relatively high at the wall according to (E11). Moreover, its gradient $\partial \sqrt { \langle \theta '^{+} \theta '^{+} \rangle } / \partial x_3$ along the normal direction to the wall is zero, corresponding to the Neumann boundary condition. The maximum r.m.s. intensity is characterized by a turbulent peak that is more pronounced for high Prandtl numbers than for low Prandtl numbers but, surprisingly, it remains roughly of the same order of order of magnitude for both cases when moving away from the wall to the central region of the channel. The turbulent peak moves towards the wall region as the Prandtl number increases from $P_r=0.01$ to 10. This section therefore confirms that the impact of the wall scalar fluctuation on the scalar field is appreciable within the vicinity of the wall but fades rapidly as the wall distance increases. For the Prandtl number unity, the distributions of the r.m.s. scalar variance in the absence of wall scalar fluctuations and the one for the turbulence energy are roughly the same because of the analogy between the velocity and scalar fields (Kim et al. 1987; Kim & Moin 1989).

Figure 4. The r.m.s. of the scalar variance $\theta ^+_{rms} = \sqrt { \langle \theta '^{+} \theta '^{+} \rangle }$ vs the wall unit distance. Here, $\theta '_{w}=0$ (case I), $\bullet$ (maroon); $q'_w=0$ (case II), $\blacktriangle$ (blue) . Panels show (a) $P_r=0.01$; (b) $P_r=0.1$; (c) $P_r=1$; (d) $P_r=10$; $R_\tau =395$.

4.5. Turbulent scalar fluxes

Figure 5 describes the profile of the turbulent scalar fluxes $\tau ^+_{i\theta }$ in the streamwise and normal directions, respectively, for the Prandtl numbers $P_r=0.01$, 0.1, 1 and 10. As expected, the scalar flux is of higher magnitude in the streamwise direction than in the normal direction to the wall. This is due to the fact that the streamwise velocity $u'_1$ is much larger than the wall-normal velocity $u'_3$ because of the wall boundary planes. The turbulent scalar fluxes increase in the immediate vicinity of the wall with the increase of the Prandtl number and return to zero when moving to the centre of the channel. Because of the change of sign of the velocity component, the streamwise scalar flux $\tau _{1\theta }$ is symmetric whereas the wall-normal scalar flux $\tau _{3\theta }$ is antisymmetric vs the wall distance. Contrary to the r.m.s. intensity of the scalar variance, which is significantly modified in the wall region, it appears that the turbulent scalar fluxes look similar from each other in appearance whatever the boundary condition applied at the wall. The reason is probably due to the fact that the transport equation for the turbulent scalar fluxes does not contain explicitly the scalar variance in their source terms and this is reinforced by the no-slip velocity boundary condition applied at the wall in all cases that imposes the correlation $\langle u'^{+}_i \theta '^{+} \rangle$ to return to zero at the wall, even if the scalar fluctuation $\theta '^{+}$ is non-zero, as verified in (E7) and (E9).

Figure 5. Turbulent scalar flux in the streamwise direction (a) $\tau ^+_{1\theta } = \left \langle u'^{+}_1 \theta '^{+} \right \rangle$ and in the normal direction, (b) $|\tau ^+_{3 \theta }| =- \left \langle u'^{+}_3 \theta '^{+} \right \rangle$ vs the wall unit distance. For $\theta '_{w}=0$ (case I): $P_r=0.01$, $\blacksquare$ (green); $P_r=0.1$, $\blacktriangle$ (blue); $P_r=1$, $\blacktriangledown$ (cyan); $P_r=10$, $\bullet$ (maroon). For $q'_w=0$ (case II): $P_r=0.01$, $\square$ (green); $P_r=0.1$, $\vartriangle$ (blue); $P_r=1$, $\triangledown$ (cyan); $P_r=10$, $\circ$ (maroon). Here, $R_\tau =395$.

4.6. Correlation coefficient

Figures 6 and 7 display the profiles of the correlation coefficients $R_{1 \theta }$ and $R_{3 \theta }$, respectively, vs the wall unit distance. The profile of the correlation coefficient $R_{1 \theta }$ is symmetric, leading to a zero gradient value in the centreline $(\partial R_{1 \theta }/\partial x_3)_c =0$ whereas $R_{3 \theta }$ is anti-symmetric, implying that it reduces to zero in the centreline $(R_{3 \theta })_c=0$. For this reason, the correlation coefficient $R_{1 \theta }$ is then of higher magnitude than $R_{3 \theta }$ almost everywhere in the channel, in accordance also with the preceding observation made in figure 5 for the turbulent fluxes showing that the scalar fluctuation is more correlated with the streamwise velocity fluctuation $u'_1$ than with the wall-normal velocity fluctuation $u'_3$. The correlation $R_{1 \theta }$ reaches its maximum value in the wall region at $P_r=1$ because of the analogy between the velocity and scalar fields. This similarity becomes, however, weaker as the Prandtl number departs from unity because the dynamic of the velocity in (2.2) and the one of the passive scalar in (2.4) are governed by different mechanisms. Equations (E16a,b) and (E17a,b) indicate that these coefficients take on finite values different from zero at the wall. The curves associated with $R_{1 \theta }$ and $R_{3 \theta }$ depart from each other at the wall but coincide away from the wall, showing once again that the effect of the wall scalar fluctuation is substantial but remains still limited to the vicinity of the wall. Both correlation coefficients for case I are of higher magnitude than their corresponding coefficients for case II. Physically, this outcome means that the passive scalar is less correlated with the velocity for the isoflux boundary condition (case II) than for the isoscalar boundary condition (case I) because the velocity fluctuation is always zero at the wall but not anymore the scalar fluctuation, the similitude between the velocity and scalar fields being lost. In case II, the scalar fluctuations are mainly diffused from the production zone (the maximum of production located somehow away from the very wall) towards the wall losing their interaction with the velocity fluctuations.

Figure 6. Correlation coefficient of the streamwise turbulent scalar flux $R_{1 \theta }$ vs the wall unit distance. Here, $\theta '_{w}=0$ (case I), $\bullet$ (maroon); $q'_w=0$ (case II), $\blacktriangle$ (blue). Panels show (a) $P_r=0.01$; (b) $P_r=0.1$; (c) $P_r=1$; (d) $P_r=10$; $R_\tau =395$.

Figure 7. Correlation coefficient of the wall-normal turbulent scalar flux $R_{3 \theta }$ vs the wall unit distance. Here, $\theta '_{w}=0$ (case I), $\bullet$ (maroon); $q'_w=0$ (case II), $\blacktriangle$ (blue). Panels show (a) $P_r=0.01$; (b) $P_r=0.1$; (c) $P_r=1$; (d) $P_r=10$; $R_\tau =395$.

4.7. Budget for the turbulent kinetic energy $k$

The budget equation for the turbulent kinetic energy $k^+ = \left \langle u'^+_j u'^+_j \right \rangle /2$ is briefly presented in figure 8 showing the profiles of the production $P^{+}_{k}$, the dissipation rate $\epsilon ^+$, the pressure diffusion $\varPhi ^+_{k}$, the turbulent $T^+_{k}$ and the molecular diffusion $d^+_{k}$. These terms appearing in the transport equation for $k$, recalled in Appendix D, are non-dimensionalized by the factor $u^4_\tau /\nu$. As known, $P^+_k$ balances with $\epsilon ^+$ away from the wall but not in the vicinity of the wall where the molecular diffusion becomes important and even equal to the dissipation $(d^+_k)_w = (\epsilon ^+)_w$ while the turbulent diffusion $T^+_k$ due to the triple velocity correlations is zero but gradually increases and decreases when moving away from the wall. The pressure diffusion $\varPhi ^+_k$ remains negligibly small everywhere in the channel compared with the other terms.

Figure 8. Budget of the transport equation for the turbulent kinetic energy $k^+= \left \langle u'^+_j u'^+_j \right \rangle /2$ vs the wall unit distance. Mesh $M_2$. $\bullet$ (maroon), $P^{+}_{k}$; $\blacktriangle$ (red), $\epsilon ^+$; $\blacklozenge$ (blue), $d^+_{k}$; $\blacksquare$ (green), $T^+_{k}$; $\blacktriangledown$ (cyan), $\varPhi ^+_{k}$. Here, $R_\tau =395$.

4.8. Budget for the scalar variance $k_\theta$

Figure 9 exhibits the budget of the transport equation for $k^+_\theta$ given in (3.1). The budget terms are non-dimensionalized by the factor $u^2_\tau \theta ^2_\tau /\nu$.

Figure 9. Budget of the transport equation for the scalar variance $k^+_{\theta } = \sqrt { \langle \theta '^{+} \theta '^{+} \rangle }/2$ vs the wall unit distance. For $\theta '_{w}=0$ (case I): $\bullet$, $P^+_{\theta }$; $\blacktriangle$ (red), $\epsilon ^+_\theta$; $\blacklozenge$ (blue), $d^+_{\theta }$; $\blacksquare$ (green), $T^+_{\theta }$. For $q'_w=0$ (case II): $\circ$, $P^+_{\theta }$; $\triangle$ (red), $\epsilon ^+_\theta$; $\lozenge$ (blue), $d^+_{\theta }$; $\square$ (green), $T^+_{\theta }$. Panels show (a) $P_r=0.01$; (b) $P_r=0.1$; (c) $P_r=1$; (d) $P_r=10$; $R_\tau =395$.

4.8.1. Case I ($\theta '_{w}=0$)

The different terms appearing in this equation are the production, where the main contribution is here $P^{+(1)}_{\theta } = - \tau ^+_{3 \theta } \partial \left \langle \theta ^+ \right \rangle / {\partial x^+_3 }$ the dissipation-rate $\epsilon ^+_\theta = \langle ({\partial \theta '^+} / {\partial x^+_j} ) ({\partial \theta '^+} / {\partial x^+_j} )\rangle / P_r$, where among the different terms, the main contribution is given by $\langle ({\partial \theta '^+} / {\partial x^+_3} ) ({\partial \theta '^+} / {\partial x^+_3} ) \rangle$, the molecular diffusion $d^+_{\theta }$ and the turbulent diffusion $T^+_{\theta }$. The peaks of high thermal production are somewhat flattened for low Prandtl numbers but become very sharp for high Prandtl numbers and also move closer to the wall as the Prandtl number increases. Overall, these terms are of higher magnitude for high Prandtl numbers. The dominant processes in the central region of the channel are always the production as a gain term and the dissipation rate as a sink term that balance each other at all Prandtl numbers. The production increases in the near-wall region and presents a peak of intensity and then decreases gradually whereas the dissipation rate reaches its maximum only at the wall. The viscous and turbulent diffusion terms are of appreciable magnitude in the very near-wall region, showing that they both play an important role in the transfer of the passive scalar. This outcome suggests therefore that the scalar conduction cannot be neglected throughout the flow field. As known, the molecular diffusion and dissipation processes are dominant in the immediate vicinity of the wall and furthermore, $d^+_{\theta }$ and $\epsilon ^+_{\theta }$ take on the same value at the wall, $(d^+_{\theta })_w= (\epsilon ^+_\theta )_w$ as analytically demonstrated in Appendix E using the Taylor series expansion in space. The turbulent diffusion becomes increasingly important as the Prandtl number increases while the viscous turbulent diffusion, on the contrary, becomes smaller. For instance, in the case of a high Prandtl number $P_r=10$, the turbulent diffusion maxima even reach a value comparable to that of the dissipation rate but it becomes almost negligible at $P_r=0.01$. As a whole, some of these results compare well with previous DNS data although the Reynolds number is here higher (Kawamura et al. 1998, 1999), and with the recent DNS performed by Alcántara-Ávila et al. (2018) and Alcántara-Ávila & Hoyas (2021) for the Reynolds number $R_\tau = 400$ and the Prandtl numbers considered here. It is of interest to remark that the budget for the half-scalar variance $k_{\theta }$ in figure 9(c) is quite similar to the budget for the turbulence kinetic energy $k$ in figure 8, suggesting that the physical mechanisms governing these equations (3.1) and (D1), respectively, are of the same nature and magnitude but only at the Prandtl number unity.

4.8.2. Comparisons between case I ($\theta '_{w}=0$) and case II ($q'_w=0$)

The curves for case I and case II, plotted on the same figures, exhibit significant changes in the immediate vicinity of the wall but they almost overlap when moving away from the wall and practically they cannot be distinguished. In particular, the dissipation-rate term $\epsilon ^+_{\theta }$ and the viscous diffusion term $d^+_{\theta }$ are still of the same order of magnitude but they are significantly smaller at the wall in case II. This outcome must be attributed to the Neumann boundary condition implying that the correlation $\langle ({\partial \theta '^+} / {\partial x^+_3} ) ({\partial \theta '^+} / {\partial x^+_3} ) \rangle$ reduces to zero at the wall. Besides the dissipation and the viscous diffusion processes that are strongly modified in the near-wall region, the production term $P^+_{\theta }$ and the turbulent diffusion term $T^+_{\theta }$ are roughly the same. These terms remain relatively unaffected by the type of boundary condition, although the turbulent diffusion $T^+_{\theta }$ is, however, slightly attenuated in the near-wall region. As these present simulations lead to new results for the Reynolds and Prandtl numbers considered, the only source of possible comparison is with the previous data of Lu & Hetsroni (1995) for $R_\tau =180$, $P_r=0.72$ and Tiselj et al. (2001) for $R_\tau =171$, $P_r=1$, showing a qualitative agreement for $P_r$ equal or close to unity. In term of turbulence modelling for RANS and LES simulations making use of transport equations for $k_\theta$ and $\epsilon _\theta$, it can be finally mentioned that the boundary conditions for $k_\theta$ and $\epsilon _\theta$ must be properly specified whether wall scalar fluctuations $\theta '_w$ are considered or not in engineering applications (Sommer, So & Zhang 1994).

4.9. Budget for the turbulent scalar fluxes $\tau _{i \theta }$

Figures 10 and 11 exhibit the budget of the transport equation for $\tau _{i \theta }$ given in (3.2). The budget terms are non-dimensionalized by the factor $u^3_\tau \theta _\tau /\nu$.

Figure 10. Budget of the transport equation for the streamwise turbulent scalar flux $\tau ^+_{1\theta } = \langle u'^{+}_1 \theta '^{+} \rangle$ vs the wall unit distance. For $\theta '_{w}=0$ (case I): $\bullet$, $P^{+}_{1 \theta }$; $\blacktriangle$ (red), $\epsilon ^+_{1 \theta }$; $\blacklozenge$ (blue), $d^+_{1 \theta }$; $\blacksquare$ (green), $T^+_{1 \theta }$; $\blacktriangledown$ (cyan), $\varPi ^+_{1 \theta }$. For $q'_w=0$ (case II): $\circ$, $P^+_{1 \theta }$; $\triangle$ (red), $\epsilon ^+_{1 \theta }$; $\lozenge$ (blue), $d^+_{1 \theta }$; $\square$ (green), $T^+_{1 \theta }$; $\triangledown$ (cyan), $\varPi ^+_{1 \theta }$. Panels show (a) $P_r=0.01$; (b) $P_r=0.1$; (c) $P_r=1$; (d) $P_r=10$; $R_\tau =395$.

Figure 11. Budget of the transport equation for the wall-normal turbulent scalar flux $\tau ^+_{3\theta } = \langle u'^{+}_3 \theta '^{+} \rangle$ vs the wall unit distance. For $\theta '_{w}=0$ (case I): $\bullet$, $P^{+}_{3 \theta }$; $\blacktriangle$ (red), $\epsilon ^+_{3 \theta }$; $\blacklozenge$ (blue), $d^+_{3 \theta }$; $\blacksquare$ (green), $T^+_{3 \theta }$; $\blacktriangledown$ (cyan), $\varPi ^+_{3 \theta }$. For $q'_w=0$ (case II): $\circ$, $P^+_{3 \theta }$; $\triangle$ (red), $\epsilon ^+_{3 \theta }$; $\lozenge$ (blue), $d^+_{3 \theta }$; $\square$ (green), $T^+_{3 \theta }$. $\triangledown$ (cyan), $\varPi ^+_{3 \theta }$. Panels show (a) $P_r=0.01$; (b) $P_r=0.1$; (c) $P_r=1$; (d) $P_r=10$; $R_\tau =395$.

4.9.1. Streamwise turbulent scalar flux $\tau _{1 \theta }$ case I ($\theta '_{w}=0$)

According to this equation, the production $P^{+}_{1 \theta }$ is decomposed into the sum of three terms as $P^{+(1)}_{1 \theta }+P^{+(2)}_{1 \theta } +P^{+(3)}_{1 \theta }$ where the main contributions are attributed to the two terms $P^{+(2)}_{1 \theta } =- \tau ^+_{1 3} \partial \langle \theta ^+ \rangle /\partial x^+_3$ and $P^{+(3)}_{1 \theta } = - \tau ^+_{3 \theta } \partial \langle u^+_1 \rangle /\partial x^+_3$ involving the interaction between the turbulent shear stress and the mean scalar gradient as well as the interaction between the wall-normal turbulent scalar flux and the mean streamwise velocity gradient. The dissipation term $\epsilon ^+_{1 \theta }$ remains large all over the channel cross-section due to the high correlation of the fluctuating velocity-scalar gradients $\langle (\partial u'^+_1 / \partial x^+_j ) \, ( \partial \theta '^+ / \partial x^+_j )\rangle$. By the way, it is found that the rough assumption of isotropic dissipation $\epsilon ^+_{1 \theta }=0$ often used in RANS closures does not hold in the near-wall region. A more appropriate order of magnitude would be $\epsilon ^+_{i\theta } = ({\lambda }/{l}) \tau ^+_{i \theta }$ as suggested by an extended development along Lumley's principles recalled in Appendix B. It appears also that the production and dissipation terms $P^{+}_{1 \theta }$ and $\epsilon ^+_{1 \theta }$ are dominant in the whole region except in the immediate vicinity of the wall where the production reduces to zero. Figure 10 reveals that the production term $P^{+}_{1 \theta }$ balances with the sum of the dissipation-rate term $\epsilon ^+_{1 \theta }$ and the scalar-pressure gradient correlation term $\varPi ^+_{1 \theta } = - \langle \theta '^+ { \partial p'^+} / {\partial x^+_1} \rangle$ away from the wall where both the molecular and turbulent diffusion terms are negligible. The molecular diffusion term $d^+_{1 \theta }$ composed by the sum of $d^{+(1)}_{1 \theta }$ and $d^{+(2)}_{1 \theta }$ becomes large only in the near-wall region for all Prandtl numbers. The turbulent diffusion term $T^+_{1 \theta }$ is very small for low Prandtl numbers but it gradually increases as the Prandtl number increases from 0.01 to 10. This result arises from the high correlation between the velocity fluctuation $u'_i$ and the scalar fluctuation $\theta '^+$ acting in the triple correlation term $\langle u'^+_1 u'^+_3 \theta '^+ \rangle$ in the same way as for the double correlation term $\langle u'^+_i \theta '^+ \rangle$ (see figure 5). As demonstrated in Appendix E, the molecular diffusion is compensated by the scalar dissipation rate at the wall $(d^+_{1 \theta } )_w= (\epsilon ^+_{1 \theta } )_w$.

4.9.2. Comparisons between case I ($\theta '_{w}=0$) and case II ($q'_w=0$)

The results returned by the DNSs with wall scalar fluctuations for the streamwise turbulent scalar flux present the same similarities as in the preceding results obtained for the scalar variance $k_\theta$. Indeed, it is found that the curves corresponding to the different terms appearing in (3.2) depart further in the wall region but start to merge as the wall distance increases from the wall. In particular, the dissipation-rate and molecular diffusion terms are considerably attenuated at the wall when wall fluctuations are present, unlike their respective counterparts computed with no wall scalar fluctuation. But the production and turbulent diffusion terms remain again almost unchanged, whatever the type of boundary condition.

4.9.3. Wall-normal turbulent scalar flux $\tau _{3 \theta }$ case I ($\theta '_{w}=0$)

The production term $P^{+}_{3 \theta }$ appearing in (3.2) is here mainly governed by the two terms $P^{+(2)}_{3 \theta } =- \tau ^+_{3 3} \partial \langle \theta ^+ \rangle /\partial x^+_3$ and $P^{+(3)}_{3 \theta } = - \tau ^+_{3 \theta } \partial \langle u^+_3 \rangle /\partial x^+_3$ involving the interaction of the wall-normal turbulent stress $\tau ^+_{33}$ and the mean scalar gradient as well as the interaction between the wall-normal turbulent scalar flux and the mean normal velocity gradient. The scalar dissipation rate $\epsilon ^+_{3 \theta }$ is caused by the correlation of the fluctuating velocity-scalar gradients $\langle ( \partial u'^+_j / \partial x^+_j ) \, ( \partial \theta '^+ / \partial x^+_j ) \rangle$ while the scalar-pressure gradient correlation term is computed as $\varPi ^+_{3 \theta } = - \langle \theta '^+ { \partial p'^+} / {\partial x^+_3} \rangle$. At all Prandtl numbers, the production is here negative, while the dissipation rate as well as the pressure gradient correlation terms are positive. It appears that the scalar-pressure gradient correlation term $\varPi ^+_{3 \theta }$ gradually increases as the Prandtl number increases so that it contributes more significantly to the budget. Contrary to the findings obtained for the streamwise turbulent scalar flux $\tau _{1 \theta }$ in § 4.9.1 showing that the budget for $\tau _{1 \theta }$ remains relatively unaffected by the Prandtl number, it is found in the present case that the budget for $\tau _{3 \theta }$ highly depends on the Prandtl number, as shown in figure 11. In particular, for $P_r=0.01$, the production term $P^{+}_{3 \theta }$ and the scalar dissipation term $\epsilon ^+_{3 \theta }$ accounting for eddies of larger scale are prominent, the correlation term $\varPi ^+_{3 \theta }$, the molecular $d^+_{3 \theta }$ including $d^{+(1)}_{3 \theta }$ and $d^{+(2)}_{3 \theta }$, as well as the turbulent diffusion $T^+_{3 \theta }$ are entirely negligible, whereas for $P_r=10$, another situation occurs, i.e. the production term $P^{+}_{3 \theta }$ balances with the sum of the scalar-pressure gradient correlation term $\varPi ^+_{3 \theta }$, the turbulent diffusion $T^+_{3 \theta }$ and the molecular diffusion $d^+_{3 \theta }$, while the scalar dissipation rate $\epsilon ^+_{3 \theta }$ is very small except, however, in the viscous sublayer near the wall region where it is appreciable. For $P_r=0.1$, the pressure-scalar gradient correlation and the dissipation terms which are the dominant sink terms balance with the production term. But for $P_r=1$, the production term mainly balances with the sum of the pressure gradient correlation and turbulent diffusion terms while the dissipation-rate and viscous diffusion terms are rather negligibly small. These results imply that the modelling of the dissipation rate $\epsilon ^+_{3 \theta }$ must be undertaken for Prandtl numbers smaller than unity and that the assumption of isotropic dissipation $\epsilon ^+_{3 \theta }=0$ is verified only for Prandtl numbers greater than unity. An approach like the one suggested in the previous section could be fulfilled. In addition, the contribution of the scalar pressure gradient correlation term $\varPi ^+_{3 \theta }$ appears essential to take into account for the budget of (3.2), at least at high Prandtl numbers, so that it cannot be ignored in conventional turbulence modelling. Hence, $\varPi ^+_{3 \theta }$, or more formally, the pressure-scalar correlation term ${\varPhi ^+_{i \theta } }$ appearing in (3.3) must be carefully modelled in the framework of second moment closure (Schiestel 2008; Hanjalic & Launder 2011). More generally, these results associated with the budget of the transport equation for $\tau _{3 \theta }$ suggest that RANS or LES modelling of (3.2) must account for the influence of Prandtl number.

4.9.4. Comparisons between case I ($\theta '_{w}=0$) and case II ($q'_w=0$)

At a first sight, no apparent difference between these curves associated with cases I and II, respectively, can be observed from the wall to the centreline of the channel. This is probably related to the fact that all these curves collapse near the wall.

4.10. Structure of the scalar fields

4.10.1. Contours of the instantaneous scalar field $\theta$

Figure 12 shows the contours plots of the instantaneous scalar field $\theta$ in the $(x_1,x_3)$ mid-plane, respectively, for the Prandtl numbers $P_r=0.01$, 0.1, 1 and 10; $\theta '_w =0$ (case I). The case $q'_w=0$ is not illustrated because the visualization of the contour lines is almost identical, at least when the Prandtl number is moderate near unity or high like 10. It requires a very enlarged view of the contour plots to shed light on some perceptible differences, the only modification been the intensity of the scalar fluctuations that is enhanced. The interesting finding, however, is that the topology of these structures considerably changes as the Prandtl number increases from $P_r=0.01$ to 10. These structures get thinner as the Prandtl number increases because of the important decrease of the Batchelor length scale $\eta _\theta$ according to the power law $\eta _\theta =\eta _\kappa /P^{1/2}_r$. These eddies are relatively smooth and organized at the lower Prandtl number $P_r=0.01$ but become more and more chaotic as the Prandtl number increases from 0.01 to 10, the gradients being sharper. This observation suggests that the scalar field is almost large scale dominant for $P_r=0.01$ but highly turbulent with fragmented structures for $P_r=10$. As the Prandtl number increases, it is possible to clearly identify substantial detachment of swirling vortex elements growing from the boundary layer to the central region of the channel along the normal direction to the wall.

Figure 12. Contours of the instantaneous scalar field $\theta$ in the $(x_1,x_3)$ mid-plane. Here, $\theta '_{w}=0$ (case I): (a) $P_r=0.01$; (b) $P_r=0.1$; (c) $P_r=1$; (d) $P_r=10$.

Figure 13 displays the contour plots of the instantaneous streamwise velocity $u_1$ in the $(x_1,x_3)$ mid-plane for illustrating the analogy with the scalar field showed in figure 12(c) at $P_r=1$. Indeed, as can be observed, a great similarity exists between the streamwise velocity field and the passive scalar field, showing almost identical eddies, even if the interfaces of $\theta$ are sharper than those of $u_1$ (Antonia, Abe & Kawamura 2009; Pirozzoli et al. 2016). In addition, it is of interest to see that an increase in Prandtl number leads to more pronounced wall-attached eddies (Marusic & Monty 2019) which is consistent with the appearance of mean logarithmic shape of $\langle \theta \rangle$ seen in figure 3. Figure 14 is a snapshot view of the scalar field $\theta$ in the $(x_1,x_2)$ plane at the wall distance $x^+_3 =10$. This figure highlights some streaky structures that are elongated in the streamwise direction for the Prandtl number $P_r=1$ in agreement with the experimental observation of Iritani, Kasagi & Hirata (1983) as well as the DNS of Li et al. (2009). They are still clearly visible and even more elongated in the streamwise direction for the Prandtl number $P_r=10$ and appear randomly in space. But on the other hand, these streaky structures are less pronounced at the Prandtl number $P_r=0.1$ and even tend to entirely disappear for the Prandtl number $P_r=0.01$. More generally, considering that the turbulent passive scalar is purely transported and diffused by the flow without any reciprocal interaction, this scalar can be viewed as a marker which provides additional information or data about the purely dynamic turbulent flow. For instance the different plots in figure 12 are dealing with the same and unique dynamic structures but the scalar allows us to capture (and picture) large structures at low Prandtl numbers and fine structures at high Prandtl numbers, a sort of ‘integral spectral’ decomposition.

Figure 13. Contours of the instantaneous streamwise velocity field $u_1$ in the $(x_1,x_3)$ mid-plane.

Figure 14. Snapshot view of the scalar field $\theta$ in the $(x_1,x_2)$ plane at the wall distance $x^+_3 =10$. Here,$\theta '_{w}=0$ (case I): (a) $P_r=0.01$; (b) $P_r=0.1$; (c) $P_r=1$; (d) $P_r=10$.

4.10.2. Contours of the instantaneous wall-normal scalar flux $\bar {\tau }_{3\theta }$

Figure 15 displays the contours of the wall-normal scalar flux ${\bar {\tau }_{3\theta } } = \overline {u'^{+}_3 \theta '^{+} }$ for the different Prandtl numbers. The positive and negative regions of $\bar {\tau }_{3 \theta }$ are characterized by a high correlation between the velocity $u'_3$ and the scalar $\theta '$ in the near-wall region for all Prandtl numbers. These contours lines describe more or less intermittent regions (Kim & Moin 1989; Abe et al. 2004).

Figure 15. Contours of the instantaneous wall-normal scalar flux $\bar {\tau }_{3\theta } = \overline {u'^{+}_3 \theta '^{+} }$ in the $(x_1,x_3)$ mid-plane. Here, $\theta '_{w}=0$ (case I): (a) $P_r=0.01$; (b) $P_r=0.1$; (c) $P_r=1$; (d) $P_r=10$.

5. Modelling of the passive scalar transport equations

5.1. The $k_\theta -\epsilon _\theta$ transport equations

In the RANS methodology, computation of turbulent flows including scalar fields needs to solve the transport equation for the turbulent energy $k$ or the transport equation for the Reynolds stress $\tau _{ij} = \left \langle u'_i u'_j \right \rangle$ in addition to the dissipation rate $\epsilon$, as well as the transport equations for the half-scalar variance $k_\theta$ and its dissipation rate $\epsilon _\theta$ that read

(5.1)$$\begin{gather} \frac{{\rm D} k_\theta } {{\rm D}t} = P_{\theta} - \epsilon_{\theta } + J_{k_{\theta}}, \end{gather}$$

(5.2)$$\begin{gather}\frac{{\rm D} \epsilon_\theta } {{\rm D}t} = c_{{\epsilon}_{\theta \theta_1}} P_{\theta} \frac{\epsilon_\theta} {k_\theta} + c_{\epsilon_{\theta k_1}} P_k \frac{\epsilon_\theta} {k} - c_{\epsilon_{\theta k_2}} \frac{\epsilon_\theta \epsilon} {k} - c_{{\epsilon}_{\theta \theta_2}} \frac{\epsilon^2_\theta} { k_\theta} + J_{\epsilon_{\theta}} , \end{gather}$$

where ${\rm D}/{\rm D}t = \partial /\partial t + \langle u_k \rangle \partial /\partial x_k$, $P_{\theta }$, $P_k$ and $J_{{k_\theta }}$, $J_{\epsilon _{\theta }}$ are the production and diffusion terms, $c_{\epsilon _{\theta \theta _1}} =1$, $c_{\epsilon _{\theta k_1}} = 1/2$, $c_{\epsilon _{\theta k_2}} = 1/2$ and $c_{\epsilon _{\theta \theta _2 }} =1.30$ (Chaouat & Schiestel 2021b). Several formulations of the transport equation for $\epsilon _\theta$ have been proposed (Newmann, Launder & Lumley 1981; Jones & Musonge 1988; Nagano & Kim 1988; Yoshizawa 1988; Shikazono & Kasagi 1996) but (5.2) is one of the most general forms obtained in the spectral space by a theoretical formalism (Chaouat & Schiestel 2021b). Because of the mathematical complexity of solving these equations arising from the occurrence of the two time scales $k/\epsilon$ and $k_\theta /\epsilon _\theta$, it is advantageous for engineering applications to compute the passive scalar to dynamic time-scale ratio $\mathcal {R}=(k_\theta \epsilon )/(k \epsilon _\theta )$ to get an estimate of the dissipation rate $\epsilon _\theta = (k_\theta \epsilon )/(\mathcal {R} k)$. This practice is usually retained for the simulation of scalar fields around $P_r$ unity considering that $\mathcal {R} \approx 1$ but the issue to address is to know whether this value can be retained or not for small and large Prandtl numbers, with and without wall scalar fluctuations. To answer this question, figure 16 describes the evolution of the ratio $\mathcal {R}$ for the Prandtl numbers $P_r=0.01$, 0.1, 1 and 10 vs the wall unit. First at all, it can be observed that this ratio is not a universal constant but is a function of both the wall unit distance and the molecular Prandtl number $\mathcal {R} = \mathcal {R} (x_3,P_r)$. For both cases, I ($\theta '_{w}=0$) and II ($q'_w=0$), it gradually decreases from high to low values in the immediate vicinity of the wall and reaches a common asymptotic value when moving away from the wall for $x^+_3 > x^+_{3m}(P_r)$ where the minimum distance $x^+_{3m}(P_r)$ is a function of the molecular Prandtl number. Moreover, its order of magnitude highly increases as the molecular Prandtl number increases from 0.01 to 10, leading to the following numerical values returned by DNS, $R(0.01) \approx 0.1$, $R(0.1) \approx 0.40$, $R(1) \approx 0.9$ and $R(10) \approx 2.5$. Because of the dissipation rate $\epsilon _\theta$ that is close to zero at the wall for the case II, $\mathcal {R}$ takes on extremely high values at $x_3=0$ in comparison with its corresponding values obtained for the case I. Physically, this observation means that the time scale of the dynamic turbulent field $k/\epsilon$ differs from the time scale of the passive scalar field $k_\theta /\epsilon _\theta$ when the molecular Prandtl number deviates from unity, suggesting once again that the analogy between these two fields is lost. In practice, interpolation between these several values can be made to compute the adequate ratio $\mathcal {R} (x_3,P_r)$ for a given molecular Prandtl number as a function also of the dimensionless distance $x^+_3$ accounting for the wall effects.

Figure 16. Dimensionless ratio $\mathcal {R}=(k_\theta \epsilon )/(k \epsilon _\theta )$ in logarithmic coordinate vs the wall unit distance. Here, $\theta '_{w}=0$ (case I), $\bullet$ (maroon); $q'_w=0$ (case II), $\blacktriangle$ (blue). Panels show (a) $P_r=0.01$; (b) $P_r=0.1$; (c) $P_r=1$; (d) $P_r=10$; $R_\tau =395$.

5.2. Turbulent fluxes $\tau _{i \theta }$

The knowledge of the turbulent scalar fluxes $\tau _{i\theta } = \left \langle u'_i \theta ' \right \rangle$ is of particular interest in mass and heat transfer. In eddy viscosity models it is computed assuming a gradient hypothesis introduced first by Daly & Harlow (1970) and often applied in practice by many authors (Hanjalic & Launder 2011):

(5.3)\begin{equation} \tau_{i\theta} ={-} \frac{\nu_t} {Pr_t } \frac{\partial \langle\theta \rangle}{\partial x_i} , \end{equation}

where $Pr_t$ is the turbulent Prandtl number and $\nu _t$ denotes the turbulent viscosity. The assessment of $\tau _{i \theta }$ must be accurate because it also stands as input in the scalar variance equation (5.1) and its dissipation-rate equation (5.2). The Prandtl number is defined itself as the ratio of the turbulent eddy viscosity $\nu _t$ to the turbulent eddy diffusivity $\sigma _t$ as

(5.4)\begin{equation} Pr_t = \frac{\nu_t} {\sigma_t} = \frac{\tau_{13}} {\tau_{3 \theta}} \frac {\partial{\langle\theta\rangle}} {\partial x_3} \left (\frac {\partial{\langle u_1\rangle}} {\partial x_3} \right )^{{-}1} . \end{equation}

Figure 17 exhibits the profile of the turbulent Prandtl number vs the wall unit distance for all DNSs. In the numbers range $P_r=0.01$ to 10, the turbulent Prandtl number reaches a decreasing asymptotic behaviour that is almost independent of the molecular Prandtl number except for $P_r=0.01$ but reveals in the near-wall region a small variation of a few per cent for the case I and a large variation for the case II. The effect of the wall scalar fluctuations is then to appreciably reduce the turbulent Prandtl number in the immediate vicinity of the wall leading to a large deviation from the asymptote, especially at low Prandtl numbers. This result validates the hypothesis of an approximately constant turbulent Prandtl number with and without wall scalar fluctuations roughly around unity away from the wall for moderate and high molecular Prandtl numbers but not for the very low Prandtl number $P_r=0.01$. From a physical point of view, the shape of these profiles is attributed to the strong conductive effects acting in low Prandtl number flows; the scalar fluctuations are then dissipated faster as the Prandtl number decreases. As the mean velocity and scalar gradient terms appearing in (5.4) are almost identical in the presence or not of the wall scalar fluctuations, the only difference between these profiles is attributed to the correlation $\tau _{3 \theta }$ appearing in the denominator of (5.4) that is modified by the wall scalar fluctuations. In practice, the channel flow is generic to many industrial applications involving heat transfer, so that the assumption of a constant turbulent Prandtl number still holds at moderate and high molecular Prandtl numbers regardless of the thermal boundary condition applied at the wall, except, however, in the near-wall region where wall profiles can be calculated to fit the DNS data. In the near-wall region, figure 17 suggests that these wall functions must gradually increase from zero to 3 ($P_r=0.01$) or roughly to unity ($P_r=0.1$, 1, 10) as a first approximation in case of wall scalar fluctuations.

Figure 17. Profiles of the turbulent Prandtl number vs the wall unit distance. Here, $\theta '_{w}=0$ (case I), $\bullet$ (maroon); $q'_w=0$ (case II), $\blacktriangle$ (blue). Panels show (a) $P_r=0.01$; (b) $P_r=0.1$; (c) $P_r=1$; (d) $P_r=10$; $R_\tau =395$.

In the framework of second moment closure (Launder et al. 1984; Schiestel 2008; Hanjalic & Launder 2011), the scalar flux $\tau _{i \theta }$ is computed by its modelled transport equation corresponding to the exact equation (3.2) that reads in a form

(5.5)\begin{equation} \frac{{{\rm D}} \tau_{i \theta} } {{\rm D}t} = P_{i \theta} + \varPhi_{i \theta} + J_{i \theta} - \epsilon_{i \theta} , \end{equation}

where $P_{i\theta }$ is the production term composed of the three terms $P^{(1)}_{i \theta }$, $P^{(2)}_{i \theta }$ and $P^{(3)}_{i \theta }$, $\varPhi _{i \theta }$ is the redistribution term associated with the pressure-scalar gradient correlation appearing in (3.3), $J_{i \theta }$ is the total diffusion term including the three terms $d^{(1)}_{i \theta }$, $d^{(2)}_{i \theta }$ and $T_{i \theta }$. In this short section, the objective is to test the basic model (5.5), some routes of modelling are also briefly suggested but without going as far as to fully validate new models, that is beyond the scope of this paper.

The Prandtl number is fixed to unity in a first step. The basic formulation of the pressure-scalar gradient term $\varPhi _{i \theta }$ used for high Reynolds number in the framework of second moment closure is (Launder et al. 1984; Schiestel 2008; Hanjalic & Launder 2011)

(5.6)\begin{equation} \varPhi_{i \theta} ={-} C^{(i)}_{1 \theta} \frac{\epsilon} {k} \tau_{i \theta} + C_{2 \theta} \tau_{j \theta} \frac{\partial \left \langle u_i \right \rangle} {\partial x_j} , \end{equation}

where $C^{(i)}_{1 \theta }$ and $C_{2 \theta }$ are constant coefficients. The modelling of the first term of (5.6) denoted $\varPhi ^{(1)}_{i \theta }$ and due to the correlation of the velocity and scalar fluctuations is inspired by the Rotta hypothesis (Monin 1965) whereas the second term $\varPhi ^{(2)}_{i \theta }$ due to the mean velocity gradient is referred to as the quasi-isotropic model, assuming that the departure from homogeneity is not too large (Launder 1976). The usual values retained for the coefficients are around $C^{(i)}_{1 \theta }=3$ ($i=1,2,3$) and $C_{2 \theta }=0.5$, although there is no precise consensus (Launder 1976, 1988; Gibson & Launder 1978; Jones & Musonge 1988). It is worth pointing out that the tensorial model closure is only generally valid if the three coefficients $C^{(i)}_{1 \theta }$ are all equal to a single true constant, otherwise the applicability is limited to the specific case studied. The previous sections have shown that $\varPhi _{i \theta }$ remains broadly unaffected by the wall scalar fluctuations and thus for its modelling (5.6). The model $\varPhi _{i \theta }$ in (5.6) is tested using the DNS data according to figures 10 and 11 that display the budgets of the transport equations for the turbulent fluxes. Figure 18(a) shows the profile of $\varPhi ^+_{i \theta }$ indicating clearly that the linear model using these constant coefficients is not at all able to match the DNS flux $\varPhi ^{+}_{1 \theta }$ that diverges near the wall region although the agreement for the DNS flux $\varPhi ^{+}_{3 \theta }$ is, however, relatively good. The disagreement with the data arises from the streamwise turbulent flux $\tau _{1 \theta }$ that exhibits a turbulent peak in the immediate vicinity of the wall while the wall-normal turbulent flux $\tau _{3 \theta }$ is very small in comparison, as shown in figure 5.

Figure 18. (a) Profile of the exact/modelled pressure-scalar gradient term $\varPhi ^+_{i \theta } = \langle p'^+ { \partial \theta '^+}/{\partial x^+_i} \rangle$ vs the wall unit distance. The DNS $\varPhi ^{+}_{1 \theta }$, $\bullet$ (blue); standard modelling (5.6), $\blacktriangledown$ (green); DNS $\varPhi ^{+}_{3 \theta }$, $\blacksquare$ (cyan); standard modelling (5.6), $\blacktriangle$ (maroon). (b) Profiles of $C^{(1)}_{1 \theta }$ $\blacktriangledown$ (blue) and $C^{(3)}_{1 \theta }$ $\blacktriangle$ (cyan) computed from (5.6) using the DNS fluxes. Here, $P_r=1$; $R_\tau =395$.

Thanks to the present simulation, it is then possible to solve (5.6) using the exact DNS terms $\varPhi ^{e}_{1 \theta }$ and $\varPhi ^{e}_{3 \theta }$ and $C_{2 \theta }=0.5$ leading to the coefficients $C^{(1)}_{1 \theta } = -k (\varPhi ^{e}_{1 \theta }-C_{2 \theta } \tau _{3\theta } \,\, \partial < u_1 >/\partial x_3 ) / (\epsilon \tau _{1\theta } )$ and $C^{(3)}_{1 \theta } = -k (\varPhi ^{e}_{3 \theta } ) / (\epsilon \tau _{3\theta } )$ that are plotted on figure 18(b). Obviously, $C^{(2)}_{1 \theta }$ is irrelevant because of transverse homogeneity. The discrepancy between these two coefficients shows clearly that the modelling (5.6) is inappropriate to simulate the turbulent channel flow subjected to constant scalar fluxes at the wall. This outcome can be compared with what was found for the modelling of the pressure strain correlation term concerning the Rotta hypothesis (Weinstock 1981, 1982). A more general formulation for the slow term $\varPhi ^{(1)}_{i \theta }$ can be then considered such as the tensorial model accounting for the quadratic term so that $\varPhi _{i \theta }$ takes the general form (Launder 1976)

(5.7)\begin{equation} \varPhi_{i \theta} ={-} \frac{\epsilon} {k} ( C^{(i)}_{1 \theta} \tau_{1 \theta} + C'^{(i)}_{1 \theta} a_{ij} \tau_{j \theta}) + C_{2 \theta} \tau_{j \theta} \frac{\partial \left \langle u_i \right \rangle} {\partial x_j} - \frac{\epsilon} {k} f_{w,\theta} C^{(i)}_{1 \theta} \tau_{k \theta} n_k n_i , \end{equation}

where $a_{ij}$ is the dimensionless anisotropic part of the Reynolds stress $a_{ij} = (\tau _{ij} - 2/3 k \delta _{ij} )/k$ and $C^{(i)}_{1 \theta }$, $C'^{(i)}_{1 \theta }$ are now some functions of several arguments of the invariant turbulent parameters or the turbulent Reynolds number (Dol, Hanjalic & Versteegh 1999), $f_{w,\theta }$ is a damping function of the near-wall model correction and $n_i$ is the wall-normal unit vector component. The function $f_{w,\theta }$ is embedded in (5.7) in the present case to satisfy the wall limiting behaviour. Setting $C'^{(i)}_{1 \theta }/C^{(i)}_{1 \theta } = \alpha$, it is then possible to compute $C^{(i)}_{1 \theta }$ from (5.7) with the DNS data of the fluxes $\varPhi ^{e}_{i \theta }$ yielding

(5.8)\begin{equation} C^{(i)}_{1 \theta} ={-} \frac{k} {\epsilon} \,\frac{\left ( \varPhi^{e}_{i \theta} - C_{2 \theta} \tau_{j \theta} \dfrac{ \partial \left\langle u_i \right\rangle} {\partial x_j} \right) } {\tau_{i \theta} + \alpha a_{ij} \tau_{j \theta} +f_{w,\theta} \tau_{k \theta} n_k n_i } . \end{equation}

The modelling of the coefficient $C^{(i)}_{1 \theta }$ is made by using the first, second and flatness invariant parameters, $A_2 = a_{ij} a_{ji}$, $A_3 = a_{ij} a_{jk} a_{ki}$ and $A=1-9(A2-A3)/8$ that has the property to reduce to zero for two-component turbulence, in particular when approaching walls, and goes to unity for isotropic stress fields because $A_2$ and $A_3$ go to zero. As a result, the model (5.7) has been calibrated with a single function of the type $C^{(i)}_{1 \theta }= 3 \,A$ (${\rm i}=1,2,3$) corresponding to the shapes of the exact functions defined in (5.8) as indicated in figure 19(b), $C_{2 \theta } = 0.5 A$, $\alpha =-0.25$ and the damping wall function $f_{w,\theta } = 5 A_2 \exp {(-A)}$. This model returns a relatively good agreement with the DNS, as shown in figure 19(a) although a small error in the magnitude is observed. This is a good compromise to get acceptable results without unnecessarily increasing the model complexity (Shikazono & Kasagi 1996). Of course, more complete testing would be necessary before recommending this choice of model coefficients in various flow applications, but this can open a pathway to further exploration. Extension to different Prandtl numbers can be undertaken by means of the use of additional functions of the molecular Prandtl number (Chaouat & Schiestel 2021b).

Figure 19. (a) Profile of the exact/modelled pressure-scalar gradient term $\varPhi ^+_{i \theta } = \langle p'^+ { \partial \theta '^+}/{\partial x^+_i} \rangle$ vs the wall unit distance. The DNS $\varPhi ^{+}_{1 \theta }$, $\bullet$ (blue); tensorial modelling (5.7), $\blacktriangledown$ (green); DNS $\varPhi ^{+}_{3 \theta }$, $\blacksquare$ (cyan); tensorial modelling (5.7), $\blacktriangle$ (maroon). (b) Profiles of $C^{(1)}_{1 \theta }$ $\blacktriangle$ (blue) and $C^{(3)}_{1 \theta }$ $\blacktriangledown$ (cyan) computed from (5.7) using the DNS fluxes. Profile of $C_{1 \theta } = 3 A$, $\bullet$ (maroon). Here, $P_r=1$; $R_\tau =395$.

The turbulent diffusion term $T_{i\theta }$ is approximated assuming a gradient law (Daly & Harlow 1970):

(5.9)\begin{equation} T_{i\theta} = \frac{\partial} {\partial x_m} \left( C_\theta \frac { k } { \epsilon } {\tau_{ml}} \frac{ \partial \tau_{i\theta} } {\partial x_l} \right) , \end{equation}

where $C_{\theta }$ is a constant coefficient set to 0.11. An extension of (5.9) is often retained as (Launder 1976)

(5.10)\begin{equation} T_{i\theta} = \frac{\partial} {\partial x_m} \left[ C_\theta \frac { k } { \epsilon } \left ( \tau_{il} \frac{ \partial \tau_{m\theta} } {\partial x_l} + \tau_{ml} \frac{ \partial \tau_{i\theta} } {\partial x_l} \right ) \right] , \end{equation}

while a more elaborate modelling that deserves interest is (Schiestel 2008; Hanjalic & Launder 2011)

(5.11)\begin{equation} T_{i\theta} = \frac{\partial} {\partial x_m} \left[ C_\theta \frac { k } { \epsilon } \left ( \tau_{il} \frac{ \partial \tau_{m\theta} } {\partial x_l} + \tau_{ml} \frac{ \partial \tau_{i\theta} } {\partial x_l} +\tau_{l \theta} \frac{ \partial \tau_{im} } {\partial x_l} \right ) \right] . \end{equation}

It is recalled that the turbulent diffusion process is not affected by the wall scalar fluctuations (see figures 10 and 11 in § 4.9). Figure 20 displays the profiles of the diffusion terms $T^{+}_{1 \theta }$ and $T^{+}_{3 \theta }$ computed using (5.9), (5.10) and (5.11), respectively, with the DNS data vs the wall unit distance. Among these modellings, it appears that (5.11) agrees best with the DNS data although the turbulent diffusion peak for $T^{+}_{1 \theta }$ is under-estimated and the one for $T^{+}_{3 \theta }$ is not at all reproduced. This analysis finally allows us to see that none of these diffusion models is really satisfactory. The molecular diffusion $d_{i \theta }$ as well as the dissipation rate $\epsilon _{i \theta }$ are generally neglected and not included in the transport equation (5.5) when the turbulence is close to an isotropy state involving the small universal scales in the Kolmogorov region of the energy spectrum. However, these models lead to an inconstant behaviour in the near-wall region where the molecular diffusion balances the dissipation rate if no wall correction is achieved and the effect of the wall scalar fluctuations (see § 4.9) must be here accounted for in the model considering that both $(d_{i \theta })_w$ and $(\epsilon _{i \theta })_w$ are close to zero. In view of this analysis, the value $(\epsilon _{i \theta })_w=0$ is suggested in (5.5) as well as $(\epsilon _{\theta })_w \approx 0$ if wall scalar fluctuations are considered. Accordingly, (5.2) could be modified to get the correct wall behaviour.

Figure 20. Diffusion term $T^{+}_{1 \theta }$: DNS, $\bullet$ (cyan); modelling (5.11), $\blacktriangle$ (maroon); modelling (5.10), $\blacklozenge$ (green); modelling (5.9), $\blacktriangleleft$ (red). Diffusion term $T^{+}_{3 \theta }$: DNS, $\blacksquare$ (cyan); modelling (5.11), $\blacktriangledown$ (maroon); modelling (5.10), $+$ (green); modelling (5.9), $\blacktriangleright$ (red). Here, $P_r=1$; $R_\tau =395$.

5.3. Turbulent fluxes $\tau _{i \theta }$ in LES

This high resolution DNS database can also be used to validate subgrid turbulence models used in LES and hybrid RANS/LES methodologies keeping in mind, however, that instead of RANS, LES usually depends on the grid spacings (Lesieur & Métais 1996; Meneveau & Katz 2000; Chaouat 2017b; Schiestel & Chaouat 2022). Strictly speaking, the RANS and LES motion equations take exactly the same mathematical form if we assume that the commutation terms arising from the filtering operation are negligible; the difference relies only on the closure of equations (Chaouat 2017a,b). Among these numerous models, the well-known Smagorinsky model with its dynamic version is the most popular model (Germano et al. 1992) to perform simulation of turbulent flows on refined grids. In the present case, focus is on the partially integrated transport model (PITM) in the framework of the second moment closure that allows us to simulate large scales of turbulent flow on relatively coarse grids (Chaouat & Schiestel 2005, 2012). This model relies on the transport equations for the subfilter turbulent stresses $\tau ^{(s)}_{ij}$ and its dissipation rate $\epsilon$ and has been recently extended to the turbulent transfer of a passive scalar including both the transport of the subfilter variance $k^{(s)}_\theta$ and its dissipation rate $\epsilon _\theta$ (Chaouat & Schiestel 2021b). These equations look like (5.1) and (5.2) where $k$ and $k_\theta$ are replaced by $k^{(s)}$ and $k^{(s)}_\theta$, respectively, but the coefficient $c_{{\epsilon }_{\theta \theta _2}}$ appearing in the destruction term of the dissipation rate $\epsilon _\theta$ is now a function $c^{(s)}_{{\epsilon }_{\theta \theta _2}}(L,\kappa _c,P_r)$ of the turbulence length scale $L=k^{3/2}/\epsilon$, the cutoff wavenumber $\kappa _c= {\rm \pi}/\varDelta$ computed using the grid size and the molecular Prandtl number. The subfilter scalar flux $\tau ^{(s)}_{i\theta }$ is modelled in analogy with (5.3) that is transposed here in LES as

(5.12)\begin{equation} \tau^{(s)}_{i\theta} ={-} \frac{c_{\tau_\theta}} {Pr_t } \tau^{(s)}_{im} \frac{k^{(s)}} {\epsilon} \frac{\partial \bar{\theta}} {\partial x_m} , \end{equation}

where $c_{\tau _\theta }=0.22$ is a numerical coefficient, $Pr_t$ is set to unity, $\bar {\theta }$ is the filtered passive scalar. Two PITM simulations of the heated channel flow (see figure 1) for case I have been performed on several meshes of coarse and medium grid resolutions if compared with DNS (see table 2), $\varDelta ^+_1=\varDelta ^+_2 \approx 60$ and $\varDelta ^+_1=\varDelta ^+_2 \approx 30$, respectively, to study the grid effect on the scalar variable. Figure 21 exhibits the profiles of the subfilter, resolved and total turbulent scalar fluxes vs the wall distance for $Pr = 1$. As a result, the sharing out of the turbulent flux between the small and resolved scales is modified according to the grid size. In the case where the grid step size decreases, a part of the scalar flux coming from the modelled zone is injected into the resolved scale, and conversely, when the grid step size increases, then a part of the energy contained in the resolved scales is removed and fed into the modelled spectral zone, but the total turbulent scalar flux remains, however, very close to DNS, confirming that the closure (5.12) is appropriate. More generally, it can be mentioned that it is therefore perfectly possible to compare directly any filtered quantity $\bar {\phi }$ of the turbulent field returned by the large eddy simulation performed on medium/coarse grids with its corresponding value $\bar {\phi }_G = G * \phi$ determined by a filtering operation as the convolution with a filter $G$ in space where $\phi$ is given by the DNS (Chaouat & Schiestel 2021a).

Figure 21. Turbulent scalar flux in the normal direction $|\tau ^+_{3 \theta }| =- \left \langle u'^{+}_3 \theta '^{+} \right \rangle$ vs the wall unit distance returned by the PITM simulation (Chaouat & Schiestel 2021b). Here, $\theta '_{w}=0$. Subfilter scale, $\blacktriangledown$ (cyan); resolved scale, $\blacktriangle$ (cyan); total scales, $\blacklozenge$ (blue); DNS, $\bullet$ (maroon). Here, $P_r=1$; $R_\tau =395$.

6. Conclusion

Thanks to the rapid increase of super-computer power, DNSs of turbulent channel flows with passive scalar transport subjected to constant scalar fluxes at the wall accounting for zero scalar fluctuation $\theta '_w=0$ (case I) and zero scalar gradient fluctuation $(\partial \theta '/ \partial x_3)_w=0$ (case II), respectively, have been performed on several meshes of high grid resolution at the Reynolds number $R_\tau =395$ and Prandtl numbers ranging from $P_r=0.01$ to 10. The evolutions of the mean scalar field $\left \langle \theta ^+ \right \rangle$, the r.m.s. scalar fluctuations $\left \langle \theta '^+ \theta '^+ \right \rangle$, the turbulent scalar flux $\tau ^+_{i \theta }$, the correlation coefficient $R_{i \theta }$ and the budgets of transport equations for $k^+_\theta$ and $\tau ^+_{i \theta }$ have been investigated in detail.

It was highlighted that the mean scalar field as well as the turbulent scalar fluxes are not affected by the wall scalar fluctuations condition in sharp contrast with the r.m.s. scalar fluctuations and the correlation coefficient $R_{i \theta }$, which are highly modified in the immediate vicinity of the wall. It has been found that the budget for $k^+_\theta$ is largely dominated by the production and the dissipation terms that balance each other at all Prandtl numbers, whereas the molecular and turbulent diffusion terms are effective in the vicinity of the wall. The dissipation-rate $\epsilon ^+_{\theta }$ and the molecular diffusion $d^+_{\theta }$ terms are highly modified in the vicinity of the wall depending on the wall scalar boundary condition. In practice, relative to case I, the very near-wall region in case II is filled with extra scalar fluctuations, the intensity of which is determined by the level of the peak of variance in the turbulence production zone. This amount of r.m.s. fluctuating scalar is blocked at the wall because of the Neumann boundary condition that cancels any diffusion flux at the wall, creating an almost uniform intensity zone. However, this extra level of scalar variance has no impact on other quantities such as $\tau ^+_{i\theta }$ because the r.m.s. itself does not appear in the source terms of the corresponding equations.

Most of these results are consistent with the findings obtained from DNS by Li et al. (2009) in a turbulent boundary layer. As expected, the Reynolds analogy is found to be satisfied when the Prandtl number is unity but only in case I because the boundary conditions for velocities and the transported scalar must be the same for the analogy to hold. As for $k^+_\theta$, the budget for $\tau ^+_{1 \theta }$ is mainly controlled by the production term $P^{+(1)}_{1 \theta }$ and the dissipation-rate term $\epsilon ^+_{1 \theta }$ away from the wall but the diffusion terms $d^+_{1 \theta }$ and $T^+_{1 \theta }$ are, however, appreciable in the vicinity of the wall. Once again, both the dissipation and molecular diffusion terms are highly influenced by the wall scalar fluctuations but the scalar-pressure gradient term $\varPi ^+_{1 \theta }$ remains unaffected by the type of boundary condition. In contrast to the budget for $\tau ^+_{1 \theta }$, the budget for $\tau ^+_{3 \theta }$ is highly dependent on the Prandtl number and is essentially governed by the production term $P^{+}_{3 \theta }$, the dissipation term $\epsilon ^+_{3 \theta }$ but also the scalar-pressure gradient correlation term $\varPi ^+_{3 \theta }$ that greatly contributes to the budget of $\tau ^+_{3 \theta }$. But surprisingly, no substantial perceptible differences are observed for the budget $\tau ^+_{3 \theta }$ with respect to the wall boundary condition.

The structures of the scalar fields have been visualized by means of the instantaneous passive scalar and wall-normal scalar flux showing that the topology of these structures is strongly modified when increasing the Prandtl number. Further arguing in this sense, the passive scalar field can be viewed as a useful marker that allows us to picture hidden properties of the dynamic field.

Finally, it has been shown how to use this DNS database to devise and calibrate scalar flux equation models in the framework of both first and second moment closures. The passive scalar to dynamic time-scale ratio $\mathcal {R}$, the turbulent Prandtl number $Pr_t$, the pressure-scalar gradient term $\varPhi _{i \theta }$, the turbulent diffusion term $T_{i\theta }$ and the dissipation $\epsilon _{i\theta }$ have been analysed and modelled from a physical point of view.

Although beyond the scope of the present study, new additional work could be pursued usefully in order to investigate the effect of the boundary conditions on the shape of spectra from a pure spectral point of view (Batchelor 1971; Tennekes & Lumley 1972; Chaouat & Schiestel 2007, 2021b). Within this framework, it should be interesting for instance to compare the velocity-scalar flux spectrum $E_{\tau _{3\theta }} (\kappa )$ with the Reynolds shear stress spectrum $E_{\tau _{13}}(\kappa )$ and examine the scale-by-scale similarity between the turbulent momentum and scalar transfers, as done by Kawata & Tsukahara (2022).

The author hopes that this high resolution DNS database will be helpful for users and researchers involved in the community of the simulation and turbulence modelling from both a theoretical and practical point of view.