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. (Reference Kawamura, Ohsaka, Abe and Yamamoto1998, Reference Kawamura, Abe and Matsuo1999), 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 Reference Luikov1968). As discussed in detail by Kasagi et al. (Reference Kasagi, Kuroda and Hirata1989), 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. Reference Kasagi, Kuroda and Hirata1989; Li et al. Reference Li, Schlatter, Brandt and Henningson2009; Chaouat & Peyret Reference Chaouat and Peyret2019). 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 Reference Errera and Chemin2013; Flageul et al. Reference Flageul, Benhamadouche, Lamballais and Laurence2015).

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 Reference Batchelor1971; Tennekes & Lumley Reference Tennekes and Lumley1972) 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 Reference Chaouat2018; Chaouat & Peyret Reference Chaouat and Peyret2019), 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 Reference Chaouat2018; Chaouat & Peyret Reference Chaouat and Peyret2019). The computational fluid dynamics (CFD) code developed by Chaouat (Reference Chaouat2011) 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.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.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. (Reference Kasagi, Kuroda and Hirata1989) for an open channel flow and Chaouat & Peyret (Reference Chaouat and Peyret2019) 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. Reference Kim, Moin and Moser1987; Kim & Moin Reference Kim and Moin1989).

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.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.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 Reference Antonia, Abe and Kawamura2009; Pirozzoli et al. Reference Pirozzoli, Bernardini and Orlandi2016). In addition, it is of interest to see that an increase in Prandtl number leads to more pronounced wall-attached eddies (Marusic & Monty Reference Marusic and Monty2019) 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 (Reference Iritani, Kasagi and Hirata1983) as well as the DNS of Li et al. (Reference Li, Schlatter, Brandt and Henningson2009). 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 Reference Kim and Moin1989; Abe et al. Reference Abe, Kawamura and Matsuo2004).

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.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 Reference Chaouat and Schiestel2021b). Several formulations of the transport equation for $\epsilon _\theta$ have been proposed (Newmann, Launder & Lumley Reference Newmann, Launder and Lumley1981; Jones & Musonge Reference Jones and Musonge1988; Nagano & Kim Reference Nagano and Kim1988; Yoshizawa Reference Yoshizawa1988; Shikazono & Kasagi Reference Shikazono and Kasagi1996) but (5.2) is one of the most general forms obtained in the spectral space by a theoretical formalism (Chaouat & Schiestel Reference Chaouat and Schiestel2021b). 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 (Reference Daly and Harlow1970) and often applied in practice by many authors (Hanjalic & Launder Reference Hanjalic and Launder2011):

(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. Reference Launder, Reynolds, Rodi, Mathieu and Jeandel1984; Schiestel Reference Schiestel2008; Hanjalic & Launder Reference Hanjalic and Launder2011), 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. Reference Launder, Reynolds, Rodi, Mathieu and Jeandel1984; Schiestel Reference Schiestel2008; Hanjalic & Launder Reference Hanjalic and Launder2011)

(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 Reference Monin1965) 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 Reference Launder1976). 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 Reference Launder1976, Reference Launder1988; Gibson & Launder Reference Gibson and Launder1978; Jones & Musonge Reference Jones and Musonge1988). 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 Reference Weinstock1981, Reference Weinstock1982). 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 Reference Launder1976)

(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 Reference Dol, Hanjalic and Versteegh1999), $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 Reference Shikazono and Kasagi1996). 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 Reference Chaouat and Schiestel2021b).

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 Reference Daly and Harlow1970):

(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 Reference Launder1976)

(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 Reference Schiestel2008; Hanjalic & Launder Reference Hanjalic and Launder2011)

(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 Reference Lesieur and Métais1996; Meneveau & Katz Reference Meneveau and Katz2000; Chaouat Reference Chaouat2017b; Schiestel & Chaouat Reference Schiestel and Chaouat2022). 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 Reference Chaouat2017a,Reference Chaouatb). Among these numerous models, the well-known Smagorinsky model with its dynamic version is the most popular model (Germano et al. Reference Germano, Piomelli, Moin and Cabot1992) 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 Reference Chaouat and Schiestel2005, Reference Chaouat and Schiestel2012). 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 Reference Chaouat and Schiestel2021b). 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 Reference Chaouat and Schiestel2021a).

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 Reference Chaouat and Schiestel2021b). 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$.