2.1. Governing equations

A DNS of the zero-pressure-gradient (ZPG) turbulent boundary layer (TBL) is performed using the pseudo-spectral code SIMSON (Chevalier et al. Reference Chevalier, Schlatter, Lundbladh and Henningson2007). The code solves the governing equations in non-dimensional form (here, written in index notation), in particular the flow and scalar variables are non-dimensionalized as

(2.1a–e)\begin{equation} \tilde{x}_{j} = \frac{x_{j}}{\delta_{0}^{*}}, \quad \tilde{U}_{j} = \frac{U_{j}}{U_{\infty}},\quad \tilde{t} = \frac{tU_{\infty}}{\delta_{0}^{*}}, \quad \tilde{P} = \frac{P}{\rho U_{\infty}^2}, \quad \tilde{\varTheta}_{i} = \frac{\varTheta_{i}-\varTheta_{i,w}}{\varTheta_{i,\infty}-\varTheta_{i,w}}, \end{equation}

where $(x_1,x_2,x_3) = (x,y,z)$ are the Cartesian coordinates in the streamwise, wall-normal and spanwise directions, respectively, and $t$ denotes the time. The length scale used for the non-dimensionalization is the displacement thickness at $x=0$ and $t=0$, denoted by $\delta _{0}^{*}$. The corresponding instantaneous velocity components are denoted by $(U_1,U_2,U_3)=(U,V,W)$ with the mean quantities identified by $(\left \langle U\right \rangle,\left \langle V\right \rangle,\left \langle W\right \rangle )$ and the fluctuations by $(u,v,w)$. Here $U_{\infty }$ is the undisturbed laminar free-stream velocity at $x=0$ and time $t=0$. The total pressure is denoted by $P$ and the density and kinematic viscosity of the fluid are represented by $\rho$ and $\nu$, respectively. In this study, four different passive scalars $(\varTheta _1,\varTheta _2,\varTheta _3,\varTheta _4)$ are simulated at different Prandtl numbers $(Pr=1,2,4,6)$, respectively. Here, $\varTheta _{i,\infty }, \varTheta _{i,w}$ correspond to the $i$th scalar concentration in the free stream and at the wall, respectively, with the mean quantities indicated by $\left \langle \varTheta _i\right \rangle$ and the corresponding fluctuations by $\theta _i$, $i\in [1,2,3,4]$. The superscript $\tilde {\cdot }$ introduced in (2.1a–e) identifies a non-dimensional variable and it shall be dropped in the non-dimensional quantities for the rest of the sections for simplicity.

The non-dimensional form of the incompressible Navier–Stokes equation and the transport equation for passive scalars are given by

(2.2)\begin{gather} \frac{\partial U_{j}}{\partial x_{j}} = 0, \end{gather}

(2.3)\begin{gather}\frac{\partial U_{j}}{\partial t} + U_{k}\frac{\partial U_{j}}{\partial x_{k}} ={-}\frac{\partial P}{\partial x_{j}} + \frac{1}{{\textit{Re}}_{\delta_0^*}}\frac{\partial^2U_{j}}{\partial x_{k}\partial x_{k}} + F_{j}, \end{gather}

(2.4)\begin{gather}\frac{\partial \varTheta_{i}}{\partial t} + U_{k}\frac{\partial \varTheta_{i}}{\partial x_{k}} = \frac{1}{{\textit{Re}}_{\delta_0^*}{\textit{Pr}}}\frac{\partial^2\varTheta_{i}}{\partial x_{k}\partial x_{k}} + F_{\varTheta_i}, \end{gather}

where ${\textit {Re}}_{\delta _0^*}$ identifies the Reynolds number based on the free-stream velocity $(U_{\infty })$ and the displacement thickness at the inlet $(\delta _0^*)$. The product of the Reynolds number based on free-stream velocity and the Prandtl number results in another non-dimensional number, called the Péclet number (${\textit {Pe}} = {\textit {Re}}_{\delta _0^*}Pr$), which measures the ratio between the scalar convective transport and the scalar molecular diffusion. Here, $F_{j}$ and $F_{\varTheta _i}$ correspond to the volume force terms for the velocity and passive scalars, respectively. The velocity–vorticity formulation of the incompressible Navier–Stokes equation is implemented in the solver as the divergence-free condition is implicitly satisfied by the formulation.

2.2. Boundary conditions

Having defined the governing equations, the problem definition is completed by providing appropriate boundary conditions. At the wall, the velocity of the fluid is the same as that of the solid surface and is given by the following no-slip and no-penetration boundary conditions:

(2.5a–c)\begin{equation} U\rvert_{y=0} = 0,\quad V\rvert_{y=0} = 0, \quad W\rvert_{y=0} = 0. \end{equation}

From the continuity equation, we also obtain

(2.6)\begin{equation} \left. \frac{\partial V}{\partial y}\right\rvert_{y=0} = 0. \end{equation}

The flow is assumed to extend to an infinite distance perpendicular to the plate, but discretizing an infinite domain is not feasible. Hence, a finite domain has to be considered, for which artificial boundary conditions have to be applied. A simple Dirichlet condition can be considered; however, the desired flow solution generally contains a disturbance that would be forced to zero. This would introduce an error due to the increased damping of the disturbances in the boundary layer (Lundbladh et al. Reference Lundbladh, Berlin, Skote, Hildings, Choi, Kim and Henningson1999). An improvement to the aforementioned boundary condition can be made by using the Neumann boundary condition given by

(2.7)\begin{equation} \left.\frac{\partial U_{j}}{\partial y}\right\rvert_{y=y_{L}} = \left.\frac{\partial \mathcal{U}_{j}}{\partial y}\right\rvert_{y=y_{L}}, \end{equation}

where $y_{L}$ is the height of the solution domain in the wall-normal direction in physical space and $\mathcal {U}_{i}$ is the laminar base flow that is chosen as the Blasius flow for the present study. For the passive scalars a constant (isothermal) wall boundary condition is applied, as given by

(2.8)\begin{equation} \varTheta_{i}\rvert_{y=0}= \varTheta_{i,w} = 0, \end{equation}

which corresponds to a vanishing thermal-activity ratio $K$. The thermal-activity ratio defines the ratio between the fluid density, thermal conductivity and specific heat capacity and the same properties of the boundary surface as defined below

(2.9)\begin{equation} K = \sqrt{\frac{\rho k C_{p}}{\rho_{w} k_{w} C_{p,w}}}. \end{equation}

Here, $\rho _{w}$, $k_{w}$ and $C_{p,w}$ correspond to the density, thermal conductivity and specific heat capacity of the wall. The isothermal wall boundary condition corresponds to the fluid that exchanges heat with the boundary surface, without modifying the wall temperature. The boundary condition in the free stream is given by

(2.10)\begin{equation} \varTheta_{i}\rvert_{y=y_{L}} = \varTheta_{i,\infty} = 1. \end{equation}

2.3. Numerical scheme

The DNS is performed with a pseudo-spectral method, where Fourier expansions are used in the streamwise and spanwise directions and Chebyshev polynomials $T_{k} (\xi )$ (on $-1\le \xi \le 1$) are used in the wall-normal direction employing the Chebyshev-tau method for faster convergence rates. The time advancement is performed using the second-order Crank–Nicholson scheme for linear terms and the third-order Runge–Kutta method for nonlinear terms, with a constant time step $\Delta t$. The maximum Courant number is set to 0.6. The nonlinear terms are calculated in physical space and the aliasing errors in the evaluation of these terms are removed by the 3/2 rule.

Since the TBL is developing in $x$, the periodic boundary condition cannot be directly used in this particular direction, which requires a specific numerical treatment. In this regard, one approach is to impose an appropriate instantaneous velocity and scalar profile at the inlet for every time step. Assuming self-similarity of the flow in the streamwise direction, Lund, Wu & Squires (Reference Lund, Wu and Squires1998) proposed a rescaling–recycling method to generate the required inlet profiles based on the solution downstream. An alternative approach is the addition of the fringe region downstream of the physical domain to retain the periodicity in the streamwise direction as described by Bertolotti, Herbert & Spalart (Reference Bertolotti, Herbert and Spalart1992); Nordström, Nordin & Henningson (Reference Nordström, Nordin and Henningson1999). In this method, the disturbances are damped, and the flow is forced from the outflow of the physical domain to the same profile as the inflow. The fringe technique is used in the present study, as the inflow conditions from a laminar profile followed by tripping produce natural instantaneous fluctuations for the velocity and the scalar fields (Araya & Castillo Reference Araya and Castillo2012). The fringe region is implemented by adding a volume force $(F_{j}; F_{\varTheta _i})$ to the momentum and scalar transport equations ((2.3) and (2.4)), respectively. The forcing term is given by

(2.11)\begin{gather} F_{j} = \lambda(x)\left(\mathcal{U}_{j} - U_{j}\right), \end{gather}

(2.12)\begin{gather}F_{\varTheta_i} = \lambda(x)\left(\breve{\varTheta}_{i} - \varTheta_{i}\right), \end{gather}

where $\lambda (x)$ is the strength of the forcing, which is non-zero only in the fringe region. The flow field at the inlet is the laminar Blasius profile $\mathcal {U}_{j}$ and for the scalar $\varTheta _i$ it is the linear variation with $y$ ranging from 0 to 1, denoted by $\breve {\varTheta }_i$.

2.4. Computational domain and numerical set-up

The laminar base flow is tripped by a random volume force strip $(\text {at }x/\delta _0^* = 10)$ to trigger the transition of the flow to a turbulent state. For this simulation, a three-dimensional cuboid is considered with length, height and width equal to $x_{L}, y_{L}, z_{L}$, respectively. The lower surface of the cuboid is considered as a flat plate with no-slip boundary conditions. The boundary layer grows in the considered computational domain with initial thickness denoted by $\delta _{0}^{*}$. In the streamwise direction, the computational domain terminates with the fringe region. The vertical extent of the computational domain includes the whole boundary layer and the domain height is chosen based on the free-stream boundary condition. In particular, the problem set-up in this work is similar to that studied by Li et al. (Reference Li, Schlatter, Brandt and Henningson2009).

The computational domain is discretized by $N_{x}, N_{y}$ and $N_{z}$ grid points in the streamwise, wall-normal and spanwise directions, respectively. The grid spacing is uniform in the streamwise and spanwise directions. For the wall-normal direction, the collocation points follow the Gauss–Lobatto distribution given by

(2.13)\begin{equation} y_{c} = \cos\left({\rm \pi}\frac{c}{N_{y}}\right) \quad c = 0,1,2,\ldots,N_{y}. \end{equation}

The computational box has dimensions $1000\delta _0^*\times 40\delta _0^*\times 50\delta _0^*$ in the streamwise, wall-normal and spanwise directions, respectively. The number of grid points in each direction is $3200\times 385\times 320$, correspondingly. Note that the smallest scale in the scalar fluctuations is inversely proportional to ${\textit {Pr}}^{1/2}$ (Tennekes & Lumley Reference Tennekes and Lumley1972) and, hence, the Batchelor length scale $\eta _{\varTheta _i}$ (which is analogous to the smallest scale in turbulent flow, the Kolmogorov scale $\eta$) is estimated as $\eta Pr^{-1/2}$ (Kozuka et al. Reference Kozuka, Seki and Kawamura2009). Similarly, the ratio of the largest to the smallest scales is proportional to $Re^{1.5}Pr^{0.5}$ at high ${\textit {Pr}}$ (Batchelor Reference Batchelor1959; Tennekes & Lumley Reference Tennekes and Lumley1972). In this study, an adequate grid resolution is adopted to resolve all the physically relevant scales. The Reynolds number based on free-stream velocity and displacement thickness at the inlet is ${\textit {Re}}_{\delta _0^*} = 450$ and the friction Reynolds number based on local friction velocity $(u_\tau )$ and boundary-layer thickness $(\delta _{99})$ is ${\textit {Re}}_\tau = 46$. At the outlet, the Reynolds number based on displacement thickness is ${\textit {Re}}_{\delta _0^*} = 1,580$ and ${\textit {Re}}_\tau = 396$.

Considering the friction velocity $u_{\tau }$ at the middle of the computational domain $(x/\delta _0^*=500\text {, corresponding to }{\textit {Re}}_\theta =794)$, the grid resolution in viscous units is $\Delta x^+=6.6$ and $\Delta z^+=3.3$ in the wall-parallel directions. In the wall-normal direction, we have an irregular distribution of collocation points, hence $\Delta y^+$ varies between $0.01$ and $3.5$. The results discussed at a particular streamwise position are inner scaled with respect to the local friction velocity, whereas, when the data along the streamwise direction are considered, the reference quantity at $x/\delta _0^* = 500$ is considered for inner scaling. This indicates that if a particular wall-parallel plane ($x$–$z$ plane) is sampled at a wall-normal distance (for e.g. $y^+=y_s$), considering the friction velocity at $x/\delta _0^* = 500$, the actual inner-scaled location varies along the streamwise direction, which is within $\pm 0.1y_s$ for the range of Reynolds numbers simulated in this work.

In this study, five different realizations of ZPG TBL are performed by introducing different trip forcings through modification of the random seed parameter, to obtain an ensemble average of the statistical quantities. All the different realizations are run for approximately 2400 time units $({\delta _{0}^{*}}/{U_{\infty }})$ after one flow through of initial transience, which corresponds to 1000 time units $({\delta _{0}^{*}}/{U_{\infty }})$. The converged statistics are obtained with the data corresponding to 12,000 time units $({\delta _{0}^{*}}/{U_{\infty }})$, equivalent to 600 eddy-turnover times.

3.1. Integral quantities and non-dimensional numbers

The shape factor $H_{12}$, which measures the ratio between the displacement thickness $\delta ^{*}$ and the momentum thickness $\theta$, is plotted in figure 1. The shape factor is lower in the turbulent region with increasing ${\textit {Re}}_{\theta }$ and agrees well with the experimental and numerical data for ${\textit {Re}}_{\theta } > 600$ where the TBL is fully developed.

Figure 1. Streamwise evolution of the shape factor $H_{12}$. (——, blue) Present DNS, ($\square$, blue) Simens et al. (Reference Simens, Jiménez, Hoyas and Mizuno2009), experimental data by ($\circ$, blue) Roach & Brierly (Reference Roach and Brierly1992), ($\vartriangle$, blue) Erm & Joubert (Reference Erm and Joubert1991), ($\triangledown$, blue) Purtell, Klebanoff & Buckley (Reference Purtell, Klebanoff and Buckley1981).

Figure 2 depicts the streamwise variation of the skin-friction coefficient $C_{f}$ for the present simulation and indicates that the obtained result is in good agreement with the turbulent skin-friction solution provided by Schoenherr (Reference Schoenherr1932), which is given by

(3.1)\begin{equation} C_{f} = 0.31\left[\ln^2{\left(2{\textit{Re}}_{\theta}\right)} + 2\ln{\left(2{\textit{Re}}_{\theta}\right)}\right]^{{-}1}. \end{equation}

The computed skin-friction coefficient is also in good agreement with the correlation proposed by Smits, Matheson & Joubert (Reference Smits, Matheson and Joubert1983), which is given as

(3.2)\begin{equation} C_{f} = 0.024 {\textit{Re}}_{\theta}^{{-}1/4}. \end{equation}

The trip location is at $x/\delta _0^*=10$, with a strong peak in the skin-friction coefficient followed by the transition to turbulence before $x/\delta _0^*=200$. The experimental data provided by Erm & Joubert (Reference Erm and Joubert1991) with wire tripping also closely correspond to the calculated turbulent skin-friction coefficient. It should be noted that Erm & Joubert (Reference Erm and Joubert1991) also repeated the experiments with different tripping devices and found the influence of tripping to persist until ${\textit {Re}}_{\theta } \approx 1500$. Due to this, Jiménez et al. (Reference Jiménez, Hoyas, Simens and Mizuno2010) found the experimental data to be scattered for ${\textit {Re}}_{\theta } <1070$ and also showed the scatter to decrease with ${\textit {Re}}_{\theta } > 1600$.

Figure 2. Variation of skin-friction coefficient along the streamwise direction. (——, blue) Present DNS, ($\cdots \cdots$, blue) theoretical laminar skin-friction solution, (- - -, blue) turbulent solution given by Schoenherr (Reference Schoenherr1932), ($\circ$, blue) experimental data by Erm & Joubert (Reference Erm and Joubert1991) with wire tripping.

The Stanton number measures the convective heat transfer into the fluid with respect to the thermal capacity of the fluid. The spatial evolution of the Stanton numbers for different passive scalars scaled with the square root of ${\textit {Pr}}$ (as plotted in figure 3) is very similar to the skin-friction profiles plotted in figure 2. There is a difference observed in the laminar Stanton number at $x=0$ with the present scaling. However, the individual Stanton-number profiles correspond well to the laminar solution (not shown in the figure) given by Kays & Crawford (Reference Kays and Crawford1993), which is

(3.3)\begin{equation} {\textit{St}} = \frac{0.332}{\sqrt{{\textit{Re}}_{x}}{\textit{Pr}}^{2/3}}, \end{equation}

and to the turbulent solution obtained from the Reynolds–Colburn analogy as given by Lienhard & John (Reference Lienhard and John2005). For $\varTheta _{1}$, due to the Reynolds analogy, the Stanton profile matches with the skin-friction profile scaled by a factor of two. For the passive scalars at higher Prandtl numbers, a generalization of the Reynolds–Colburn analogy can be obtained, as reported in the study by Lienhard (Reference Lienhard2020)

(3.4)\begin{equation} {\textit{St}} = \frac{Nu}{{\textit{Re}}_{x}{\textit{Pr}}} = \frac{C_{f}/2}{a_{1}+a_{2}\left({\textit{Pr}}^{a_{3}}-1\right)\sqrt{C_{f}/2}}, \end{equation}

with the values of $a_{1}=1, a_{2} = 12.8$ and $a_{3}=0.68$ provided in Lienhard & John (Reference Lienhard and John2005) and $C_f$ being expressed as

(3.5)\begin{equation} C_{f} = \frac{0.455}{\ln^2{\left(0.06{\textit{Re}}_{x}\right)}}. \end{equation}

The Stanton-number plot for the passive scalar at ${\textit {Pr}} = 6$ agrees well with the turbulent solution obtained from the Reynolds–Colburn analogy, as shown in figure 3. On the other hand, reference curves for Kays & Crawford and the Reynolds–Colburn analogy are only reported for $\varTheta _{4}$ for clarity.

Figure 3. Variation of Stanton number along the streamwise direction for different passive scalars. (——, orange) $\varTheta _{1}\,$, (——, green) $\varTheta _{2}$, (——, red) $\varTheta _{3}$, (——, blue) $\varTheta _{4}$, ($\cdots \cdots$, blue) Kays & Crawford (Reference Kays and Crawford1993) correlation corresponding to $\varTheta _4$ and (- - -, blue) Reynolds–Colburn analogy as given by Lienhard & John (Reference Lienhard and John2005) plotted for $\varTheta _{4}$, (— $\cdot$, blue) interpolation developed by Hollingsworth (Reference Hollingsworth1989) plotted for $\varTheta _{4}$.

Using the data for water with ${\textit {Pr}} = 5.9$, Hollingsworth (Reference Hollingsworth1989) developed an interpolation for Prandtl numbers from 0.7 to 5.9 assuming the critical thickness of the sub-layer to be a simple function of Prandtl number. The empirical expression is given by

(3.6)\begin{equation} {\textit{St}} = 0.02426 {\textit{Pr}}^{{-}0.895} {\textit{Re}}_{x}^{{-}0.1879 {\textit{Pr}}^{{-}0.18}}, \end{equation}

which is plotted for the scalar at ${\textit {Pr}} = 6$ in figure 3. We find that the relationship proposed by Hollingsworth (Reference Hollingsworth1989) underpredicts the Stanton number at ${\textit {Pr}} = 6$. However, the interpolation relation provides a good match with the calculated Stanton number for $\varTheta _{1}$, better than the Reynolds–Colburn analogy, which is not indicated in the plot for clarity.

Using the data obtained in the present simulation, a correlation between the Nusselt, Prandtl and Reynolds numbers can be obtained as

(3.7)\begin{equation} {\textit{Nu}}_{x} = 0.02 {\textit{Re}}_{x}^{0.828} {\textit{Pr}}^{0.514}, \end{equation}

which yields $R^{2} = 0.9985$. This is similar to the correlation proposed by Kays & Crawford (Reference Kays and Crawford1993) but for fully developed profiles in circular tubes and computed for gases

(3.8)\begin{equation} {\textit{Nu}} = 0.021 Re^{0.8} {\textit{Pr}}^{0.5}. \end{equation}

3.2. Mean velocity and scalar profiles

The mean velocity profile obtained at the streamwise location corresponding to ${\textit {Re}}_{\theta } = 670$ is shown in figure 4. The streamwise velocity profile is compared with the DNS data from Spalart (Reference Spalart1988) at ${\textit {Re}}_{\theta } = 670$ and Komminaho & Skote (Reference Komminaho and Skote2002) at ${\textit {Re}}_{\theta } = 666$. The comparison of the present data with the existing DNS results shows a good agreement in the inner region. There is a slight deviation of the mean velocity profile reported by Spalart (Reference Spalart1988) in the wake region with respect to the present data but it agrees well with the data provided by Komminaho & Skote (Reference Komminaho and Skote2002).

Figure 4. Inner-scaled mean streamwise velocity profile at ${\textit {Re}}_{\theta } = 670$. (——, blue) Present DNS, ($\circ$, blue) Spalart (Reference Spalart1988), ($\square$, blue) Komminaho & Skote (Reference Komminaho and Skote2002) at ${\textit {Re}}_{\theta } = 666$.

The mean profiles of the various passive scalars $\varTheta _{i}$ are normalized with the respective Prandtl numbers and represented in inner scaling by dividing with the friction scalar $\varTheta _{i,\tau }$ defined as

(3.9)\begin{equation} \varTheta_{i,\tau} = \frac{q_{i,w}}{\rho C_{p} u_{\tau}}, \end{equation}

where $C_{p}$ is the heat capacity of the fluid and $q_{i,w}$ is the rate of heat transfer from the wall to the fluid and is defined by

(3.10)\begin{equation} q_{i,w} = \left.-k\frac{{\rm d}\left\langle \varTheta_{i}\right\rangle}{{\rm d} y}\right\rvert_{y=0}, \end{equation}

where $k$ is the thermal conductivity of the fluid. The normalized inner-scaled mean scalar profiles are plotted at ${\textit {Re}}_{\theta } = 1070$, corresponding to ${\textit {Re}}_{\tau } = 395$, as shown in figure 5. The mean scalar profiles follow the conductive sub-layer relation ($\varTheta _{i}^+={\textit {Pr}} y^+$) and is clearly identified in the plots for $y^+ < 5$. The profiles of the passive scalars at ${\textit {Pr}}=1,2$ are compared against the channel DNS data provided by Kozuka et al. (Reference Kozuka, Seki and Kawamura2009) at the same Prandtl numbers. The comparison shows a good agreement in the near-wall and overlap regions. It should be noted that the boundary condition used by Kozuka et al. (Reference Kozuka, Seki and Kawamura2009) is the uniform-heat-flux condition as opposed to the uniform scalar boundary condition applied in this study. Note that the study by Kawamura, Abe & Shingai (Reference Kawamura, Abe and Shingai2000) showed that the Kármán constant $k_{\varTheta }$ as appearing in (3.11)

(3.11)\begin{equation} \frac{\left\langle \varTheta_{i}\right\rangle - \left\langle \varTheta_{i,w}\right\rangle}{\varTheta_{i,\tau}} = \frac{1}{k_{\varTheta}} \log y^+{+} A_{\varTheta_i}\left({\textit{Pr}}\right) , \end{equation}

varies between the Dirichlet and Neumann boundary conditions based on low-Reynolds-number simulations. Here, $A_{\varTheta _i}$ denotes the additive constant for scalar $\varTheta _i$. However, based on higher-Reynolds-number simulations, Pirozzoli, Bernardini & Orlandi (Reference Pirozzoli, Bernardini and Orlandi2016) reported that the difference in boundary condition affects the mean passive scalar profiles only in small magnitudes in the logarithmic region, although the effect is evident in the scalar fluctuation profiles reported later.

Figure 5. Normalized mean profiles of the passive scalars at ${\textit {Re}}_{\theta } = 1070$ corresponding to ${\textit {Re}}_{\tau } = 395$. (——, blue) ${\textit {Pr}}=1$, (——, orange) ${\textit {Pr}}=2$, (——, green) ${\textit {Pr}}=4$, (——, red) ${\textit {Pr}}=6$, channel DNS data by Kozuka et al. (Reference Kozuka, Seki and Kawamura2009) at ${\textit {Re}}_{\tau } = 395$ and ($\circ$, blue) ${\textit {Pr}}=1$, ($\circ$, orange) ${\textit {Pr}}=2$, ($\circ$, red) ${\textit {Pr}}=7$, ($\circ$, green) channel DNS data by Alcántara-Ávila, Hoyas & Pérez-Quiles (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2018) at ${\textit {Re}}_{\tau } = 500$ and ${\textit {Pr}} = 4$, with (- - -, grey) corresponding fit by Pirozzoli (Reference Pirozzoli2023) and (inset) ($\cdots \cdots$, black) conductive sub-layer relation, $\varTheta _{i}^+=Pr y^+$.

The scalar profile at ${\textit {Pr}} = 4$ is compared against the channel DNS data reported by Alcántara-Ávila et al. (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2018). The channel DNS data reported by Kozuka et al. (Reference Kozuka, Seki and Kawamura2009) at ${\textit {Pr}} = 7$ are used for comparison of the passive scalar at ${\textit {Pr}} = 6$, since there are no simulations reported in the literature at exactly the same Prandtl number. This gives us the opportunity to highlight the difference between the profiles at these high Prandtl numbers. Due to the difference in the considered Prandtl numbers, we observe a discrepancy in the mean velocity profile for $y^+ > 40$ in the overlap region. Since the channel data are used for the comparison, there is a difference observed near the wake region for all the cases. However, a good agreement of the profiles is observed for the inner region. Based on experimental data, semi-empirical fits were provided by Kader (Reference Kader1981) for a boundary layer with constant heat flux. In this study it was assumed that the overlap layer exhibited logarithmic variation as given in (3.11), and an empirical relation was provided to determine the additive constant $B_{\varTheta _i}$. The comparison of the mean scalar profiles against the relationships provided by Kader (Reference Kader1981) shows that the value at the wake is slightly overestimated with respect to the DNS data. The deviation for scalar $\varTheta _4$ corresponding to ${\textit {Pr}} = 6$ is approximately 3 % for $y^+ \in [100,500]$. One possible reason for this small deviation could also be the constant scalar boundary condition imposed in our simulations as opposed to the constant-heat-flux boundary conditions considered by Kader (Reference Kader1981). The comparison against Kader (Reference Kader1981) is not presented, instead, the mean scalar profiles are evaluated against the recent correlations proposed by Pirozzoli (Reference Pirozzoli2023). These correlations were developed by fitting the functional form of eddy thermal diffusivity to DNS data, resulting in predictive formulae for mean scalar profiles. Our comparison demonstrates good agreement across the different Prandtl numbers. For the passive scalar $\varTheta _4$, the maximum observed deviation is approximately 2 %.

In the overlap region, the mean scalar profiles are nearly logarithmic and are characterized as

(3.12)\begin{equation} \left\langle \varTheta_i\right\rangle^+_{\infty} - \left\langle \varTheta_i\right\rangle^+{=} -\frac{1}{k_{\varTheta}} \log \left(y/\delta^*\right) + B_{\varTheta_i}\left({\textit{Pr}}\right). \end{equation}

The mean scalar and velocity profiles are plotted in defect form in figure 6. We observe that the slope of the scalar profiles in defect form is nearly parallel indicating a constant Kármán constant for the scalar fields. The Kármán constant for the velocity is observed to be $0.41$, consistent with the values reported in Spalart (Reference Spalart1988). The diagnostic function for scalar is defined as

(3.13)\begin{equation} \varXi_\varTheta = y^+ \frac{\partial \left\langle \varTheta_i\right\rangle^+ }{\partial y^+} = \frac{1}{k_{\varTheta}}, \end{equation}

is plotted in figure 7. The Kármán constant for a scalar field is obtained as $k_{\varTheta } = 0.4$, which is in between $k_{\varTheta } = 0.33$ (Wikström Reference Wikström1998) and $k_{\varTheta } = 0.47$ (Kader Reference Kader1981). The obtained $k_{\varTheta }$ from the present simulation is the same as observed by Kawamura, Abe & Matsuo (Reference Kawamura, Abe and Matsuo1999) and close to the value of $0.41$ as reported by Li et al. (Reference Li, Schlatter, Brandt and Henningson2009). In the outer region, we observe a collapse of the profiles.

Figure 6. Mean profiles in defect form at ${\textit {Re}}_{\theta } = 1070$ corresponding to ${\textit {Re}}_{\tau } = 395$. (——, blue) ${\textit {Pr}}=1$, (——, orange) ${\textit {Pr}}=2$, (——, green) ${\textit {Pr}}=4$, (——, red) ${\textit {Pr}}=6$ and (- - -, black) velocity defect $\left \langle U_\infty \right \rangle ^+ - \left \langle U\right \rangle ^+$.

Figure 7. Variation of the diagnostic function for scalar $\varXi _\varTheta$ along the wall-normal direction at ${\textit {Re}}_{\theta } = 1070$ corresponding to ${\textit {Re}}_{\tau } = 395$. (——, blue) ${\textit {Pr}}=1$, (——, orange) ${\textit {Pr}}=2$, (——, green) ${\textit {Pr}}=4$, (——, red) ${\textit {Pr}}=6$ and (— ${\cdot }$, black) $\varXi _\varTheta =0.25;\;k_\varTheta =0.4$.

3.3. Velocity and scalar fluctuations

As shown in figure 8, the present velocity-fluctuation root-mean-squared (r.m.s.) data show a trend similar to that of the results by Jiménez et al. (Reference Jiménez, Hoyas, Simens and Mizuno2010). The r.m.s. profiles of the three velocity components are in good agreement in the inner region, while a minor difference can be observed in the outer region. Nonetheless, the peaks of the velocity fluctuations match in both position and magnitude. There is a slight offset in the plots of $p_{rms}$, which can be attributed to the small difference in the considered ${\textit {Re}}_{\theta }$ for comparison. The r.m.s. of the velocity components are also compared with the channel DNS data provided by Abe et al. (Reference Abe, Kawamura and Matsuo2004). The streamwise r.m.s. agrees well with the present DNS results and the near-wall peak value coincides with the present observations. As expected, there is a difference observed in the outer region of flow, since channel and boundary-layer flows are fundamentally different farther from the wall. Additionally, the r.m.s. of the pressure fluctuations observed in the boundary layer is different compared with the channel flow. From the present DNS data, the peak of the streamwise velocity fluctuation is found at $y^+ = 14$, corresponding in outer units to $y/\delta _{99} = 0.035$.

Figure 8. The r.m.s. velocity fluctuations and r.m.s. pressure fluctuation at ${\textit {Re}}_{\theta } = 1070$. (——, blue) Present DNS, ($\circ$, blue) Jiménez et al. (Reference Jiménez, Hoyas, Simens and Mizuno2010) at ${\textit {Re}}_{\theta } = 1100$, ($\square$, blue) channel DNS data of Abe, Kawamura & Matsuo (Reference Abe, Kawamura and Matsuo2004) at ${\textit {Re}}_{\tau }=395$.

The r.m.s. of the scalars at different Prandtl numbers are plotted in figure 9. The scalar r.m.s. profile at ${\textit {Pr}} = 1$ is similar to the streamwise velocity r.m.s. and has a higher (roughly 5 %) near-wall peak comparatively, as expected. The comparison of scalar-fluctuation profiles at ${\textit {Pr}} = 1,2$ with the channel DNS data from Kozuka et al. (Reference Kozuka, Seki and Kawamura2009) shows a good agreement in the inner and logarithmic region in addition to a good match of the peak value and wall-normal location. Despite the small difference in ${\textit {Re}}_{\tau }$, the profiles at ${\textit {Pr}} = 4$ obtained by Alcántara-Ávila et al. (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2018) show a reasonably good agreement with the present results. With increasing Prandtl number, the peak value of the scalar-fluctuation r.m.s. increases and is located closer to the wall. The scalar fluctuations decay to zero at the wall due to the isothermal boundary condition and they also decay to zero outside the boundary layer due to the absence of disturbances in the free stream.

Figure 9. Scalar fluctuation r.m.s. at ${\textit {Re}}_{\theta } = 1070$. Present DNS data for (——, blue) ${\textit {Pr}} = 1$, (——, orange) ${\textit {Pr}} = 2$, (——, green) ${\textit {Pr}} = 4$, (——, red) ${\textit {Pr}} = 6$, channel DNS data by Kozuka et al. (Reference Kozuka, Seki and Kawamura2009) at ${\textit {Re}}_{\tau } = 395$ and ($\circ$, blue) ${\textit {Pr}}=1$, ($\circ$, orange) ${\textit {Pr}}=2$, ($\circ$, red) ${\textit {Pr}}=7$, ($\circ$, green) channel DNS data by Alcántara-Ávila et al. (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2018) at ${\textit {Re}}_{\tau } = 500$ and ${\textit {Pr}} = 4$.

The obtained scalar-fluctuation r.m.s. profiles are scaled with the respective Prandtl numbers and plotted in figure 10. We observe that the lines of $\theta ^+_{i\,rms}$ for different scalars at different Reynolds numbers are parallel and not coinciding. Similar observations were also made by Alcántara-Ávila, Hoyas & Pérez-Quiles (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021) for a particular scalar at ${\textit {Pr}} = 0.71$ and they attributed the differences in the viscous-diffusion term at the wall to the increase in the slope of $\theta ^+_{i\,rms}$. The slope of $\theta ^+_{i\,rms}$ changes with Reynolds number because the peak of $\theta ^+_{i\,rms}$ varies as shown in figure 11 and the location of such a peak moves farther from the wall as the Reynolds number increases. On the other hand, the location of $\theta ^+_{i\,rms, max}$ moves closer to the wall with respect to the Prandtl number, as discussed earlier. A correlation for $\theta _{i\,rms,\, max}^+$ was obtained as in Alcántara-Ávila et al. (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021) for each passive scalar considered in the present study, as shown in figure 11. It should be noted that if the maximum value is obtained directly at the collocation point, some variations are observed along ${\textit {Re}}_{\tau }$ due to the grid resolution. In order to minimize the high-frequency variation along ${\textit {Re}}_{\tau }$, as plotted in figures 11 and 12, a simple convolution operation is performed which does not alter the obtained empirical fits. Clearly, for higher Prandtl numbers of $4,6$, the peak of the scalar-fluctuation r.m.s. tends to decrease with increasing ${\textit {Re}}_{\tau }$. When the present DNS data with ${\textit {Re}}_{\theta } \in [470,1070]$ and ${\textit {Pr}} = 1,2,4,6$ are used to find an overall variation of the peak in the scalar-fluctuation r.m.s. we find

(3.14)\begin{equation} \theta_{i\,rms, max}^+{=} 2.969 {\textit{Re}}_{\tau}^{{-}0.00858} {\textit{Pr}}^{0.571}, \end{equation}

with $R^2 = 0.99$, where $R^2$ corresponds to the coefficient of determination. It is defined as $R^2 = 1-({\varSigma (x_t - x_p)^2})/({x_{t\, {rms}}^2})$, with $x_p$ corresponding to the fitted value of $x_t$. The observed correlation shows a weak dependence on ${\textit {Re}}$ compared with ${\textit {Pr}}$. Further, from (3.14), a decaying trend of $\theta _{i\, rms, max}^+$ with ${\textit {Re}}$ is obtained, although such a trend has not been observed in the literature except for $\varTheta _3, \varTheta _4$ in the present study. The studies by Pirozzoli et al. (Reference Pirozzoli, Bernardini and Orlandi2016) and Alcántara-Ávila et al. (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021) have reported the increasing trend of $\theta _{i\,rms,max}^+$ with respect to ${\textit {Re}}$ but for a lower ${\textit {Pr}}$ than in the present work. Pirozzoli et al. (Reference Pirozzoli, Bernardini and Orlandi2016) have suggested that the attached-eddy arguments support the increase of inner peak of streamwise velocity r.m.s. with respect to ${\textit {Re}}_\tau$ due to the effect of overlying attached eddies and they assumed that the same argument applies to the passive scalars. On the other hand, we find that the inner peak of the passive scalar r.m.s. at high ${\textit {Pr}}$ does not follow the previous argumentation. It should also be noted that the range of ${\textit {Re}}_\tau$ in our present simulation is narrow compared with those of the works by Pirozzoli et al. (Reference Pirozzoli, Bernardini and Orlandi2016) and Alcántara-Ávila et al. (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021) and a more detailed investigation of this topic is necessary to make any conclusive statements.

Figure 10. Scalar fluctuations scaled by Prandtl number at (——, blue) ${\textit {Pr}} = 1$, (——, orange) ${\textit {Pr}} = 2$, (——, green) ${\textit {Pr}} = 4$, (——, red) ${\textit {Pr}} = 6$ and (——, black) ${\textit {Re}}_{\theta } = 420$, (- - -, black) ${\textit {Re}}_{\theta } = 628$, (— ${\cdot }$, black) ${\textit {Re}}_{\theta } = 830$, ($\cdots \cdots$, black) ${\textit {Re}}_{\theta } = 1070$ at the corresponding Prandtl numbers.

Figure 11. Variation of the peak of scalar fluctuation r.m.s. with respect to ${\textit {Re}}_{\tau }$, for (——, blue) ${\textit {Pr}} = 1$, (——, orange) ${\textit {Pr}} = 2$, (——, green) ${\textit {Pr}} = 4$, (——, red) ${\textit {Pr}} = 6$ including the corresponding empirical fits ($\cdots \cdots$, blue) $0.116 \log {({\textit {Re}}_{\tau })} + 2.210$ with $R^{2} = 0.976$, ($\cdots \cdots$, orange) $0.037 \log {({\textit {Re}}_{\tau })} + 3.963$ with $R^{2} = 0.620$, ($\cdots \cdots$, green) $-0.068 \log {({\textit {Re}}_{\tau })} + 6.583$ with $R^{2} = 0.726$, ($\cdots \cdots$, red) $-0.151 \log {({\textit {Re}}_{\tau })} + 8.753$ with $R^{2} = 0.894$.

Figure 12. Variation of the peak of streamwise heat flux with respect to ${\textit {Re}}_{\tau }$, for (——, blue) ${\textit {Pr}} = 1$, (——, orange) ${\textit {Pr}} = 2$, (——, green) ${\textit {Pr}} = 4$, (——, red) ${\textit {Pr}} = 6$ including the corresponding empirical fits ($\cdots \cdots$, blue) $0.250 \log {({\textit {Re}}_{\tau })} + 5.732$ with $R^{2} = 0.808$, ($\cdots \cdots$, orange) $0.064 \log {({\textit {Re}}_{\tau })} + 9.769$ with $R^{2} = 0.119$, ($\cdots \cdots$, green) $-0.1626 \log {({\textit {Re}}_{\tau })} + 14.223$ with $R^{2} = 0.328$, ($\cdots \cdots$, red) $-0.293 \log {({\textit {Re}}_{\tau })} + 16.773$ with $R^{2} = 0.513$.

The above procedure was also performed for the streamwise heat flux, as shown in figure 12. A similar behaviour was observed with the overall variation in streamwise heat flux given by

(3.15)\begin{equation} \left\langle u\theta_i\right\rangle^+_{max} = 7.74 {\textit{Re}}_{\tau}^{{-}0.006} {\textit{Pr}}^{0.401}, \end{equation}

with $R^{2} = 0.9913$.

3.4. Turbulent Prandtl number

An important parameter for scalar transport is the turbulent Prandtl number ${\textit {Pr}}_{t}$, which is defined as the ratio between turbulent eddy viscosity and turbulent eddy diffusivity

(3.16)\begin{equation} {\textit{Pr}}_{t} = \frac{\nu_{t}}{\alpha_{t}}. \end{equation}

The eddy viscosity and the eddy diffusivity arise from the Boussinesq hypothesis (Boussinesq Reference Boussinesq1877) for modelling turbulent stresses and the heat-flux vector, respectively. For parallel flows (the ones in which the velocity profile does not vary in the streamwise direction), the turbulent eddy viscosity and the turbulent eddy diffusivity are used to describe the turbulent momentum transfer and heat transfer with respect to the mean-flow conditions, in particular the mean velocity strain and temperature gradients, respectively. They are defined as

(3.17)\begin{equation} \nu_{t} ={-}\frac{\left\langle uv\right\rangle}{{\partial \left\langle U\right\rangle}/{\partial y}}, \end{equation}

and

(3.18)\begin{equation} \alpha_{t} ={-}\frac{\left\langle v\theta_i\right\rangle}{{\partial \left\langle \varTheta_i\right\rangle}/{\partial y}}. \end{equation}

It should be noted that the eddy viscosity or diffusivity does not represent a physical property of the fluid, like the molecular viscosity, but rather a property of the flow. The Reynolds analogy introduces the similarity between the turbulent momentum exchange and turbulent heat transfer in a fluid. Reynolds (Reference Reynolds1874) noted that, for a fully turbulent field, both the momentum and heat are transferred due to the motion of turbulent eddies. This yields a simpler model for the turbulent Prandtl number, where the turbulent eddy viscosity for the momentum exchange and turbulent eddy diffusivity for the scalar transport are equal, such that ${\textit {Pr}}_{t} = 1$.

Substitution of the eddy viscosity and diffusivity into the definition of turbulent Prandtl number results in

(3.19)\begin{equation} {\textit{Pr}}_{t} = \frac{\left\langle uv\right\rangle}{\left\langle v\theta_i\right\rangle}\frac{{\partial \left\langle \varTheta_i\right\rangle}/{\partial y}}{{\partial \left\langle U\right\rangle}/{\partial y}}. \end{equation}

From this definition, evaluating the turbulent Prandtl number at any point in the boundary layer would require the turbulent shear stress, turbulent heat transfer, velocity gradient and temperature gradient. Experimental investigations have limited accuracy in the simultaneous measurement of the Reynolds shear stress and wall-normal turbulent heat flux, in particular close to the wall (Araya & Castillo Reference Araya and Castillo2012). For this reason, experimental investigations like Blom (Reference Blom1970), Antonia (Reference Antonia1980) and Kays & Crawford (Reference Kays and Crawford1993) exhibit significant scatter in the data.

The variation of turbulent Prandtl number with the wall-normal distance in inner units at a given ${\textit {Re}}_{\theta } = 1070$ is reported in figure 13. We observe that the turbulent Prandtl number varies at the wall and increases with respect to molecular Prandtl number of the scalar. From figure 14, we also see the turbulent Prandtl number decays for $y^+ > 15$ and this decay becomes steeper as ${\textit {Re}}$ decreases. The turbulent Prandtl number was assumed to approach a constant value close to the wall independent of the molecular Prandtl number (Li Reference Li2007). Studies such as the ones by Kestin & Richardson (Reference Kestin and Richardson1963) and Blom (Reference Blom1970) have analysed experimentally the turbulent Prandtl number in order to assess the validity of this assumption. Such experimental investigations have reported a constant value of turbulent Prandtl number as the wall is approached. However, different experimental campaigns provided data that could not show a conclusive interpretation of the behaviour of turbulent Prandtl number close to the wall. Nonetheless, many authors have proposed different correlations based on the experimental data to predict the heat-transfer coefficient through the prediction of turbulent Prandtl number with respect to molecular Prandtl number.

Figure 13. Variation of turbulent Prandtl number with respect to wall-normal distance in inner scale at ${\textit {Re}}_{\theta } = 1070$ for passive scalar at (——, blue) ${\textit {Pr}} = 1$, (——, orange) ${\textit {Pr}} = 2$, (——, green) ${\textit {Pr}} = 4$, (——, red) ${\textit {Pr}} = 6$.

Figure 14. Variation of turbulent Prandtl number with respect to wall-normal distance in inner scale at (——, blue) ${\textit {Re}}_{\theta } = 396$, (——, orange) ${\textit {Re}}_{\theta }=420$, (——, green) ${\textit {Re}}_{\theta }=628$, (——, red) ${\textit {Re}}_{\theta }=830$, (——, purple) ${\textit {Re}}_{\theta }=1070$ for a passive scalar at ${\textit {Pr}} = 1$.

For the different passive scalars considered in the present study, the turbulent Prandtl number approaches a constant value $>1$ close to the wall in the viscous sub-layer. The plots of ${\textit {Pr}}_{t}$ exhibit a significant difference closer to the wall with respect to various molecular Prandtl numbers. There is a slight decrease in ${\textit {Pr}}_{t}$ up to $y^+ \approx 20$ and then the increase is maintained farther from the wall up to $y^+ \approx 50$, the point after which the trend steadily decreases. It is pointed out in Kays (Reference Kays1994) that the peak between $y^+ \approx 20$ and $100$ is not observed in the experimental data, an observation which was attributed to the high Reynolds number of the experiments, where DNS data are not available for comparison.

Following the discussions about a constant ${\textit {Pr}}_{t}$ in the logarithmic region, Kays & Crawford (Reference Kays and Crawford1993) proposed a correlation for the turbulent Prandtl number that is applicable to air. In this correlation, the value of ${\textit {Pr}}_{t}$ approaches a constant value of 0.85 in the logarithmic region. In the studies by Hollingsworth (Reference Hollingsworth1989), a correlation was proposed based on the measurement of the temperature profile of water at ${\textit {Pr}} \approx 6$. Again, the value of ${\textit {Pr}}_{t}$ approaches a value of 0.85 as $y^+$ is increased. This observation of constant ${\textit {Pr}}_{t}$ in the logarithmic region is not clearly observed with the present DNS data. This can be due to the low-Reynolds-number range considered in the present study. Based on the correlation proposed by Hollingsworth (Reference Hollingsworth1989), Kays (Reference Kays1994) suggested a constant value of ${\textit {Pr}}_{t} = 1.07$ for $0< y^+<5$. Indeed, if we consider only the passive scalars at ${\textit {Pr}} = 1,2$, the turbulent Prandtl number approaches a constant value of 1.07 closer to the wall. It appears that the turbulent Prandtl number close to the wall is independent of the molecular Prandtl number, as shown in the studies by Li et al. (Reference Li, Schlatter, Brandt and Henningson2009). This independence of the turbulent Prandtl number at the wall with respect to the molecular Prandtl number has also been reported in many other studies like Kong et al. (Reference Kong, Choi and Lee2000) and Jacobs & Durbin (Reference Jacobs and Durbin2000) for TBL flow, as well as in Kim & Moin (Reference Kim and Moin1989) and Kasagi, Tomita & Kuroda (Reference Kasagi, Tomita and Kuroda1992) for turbulent channel flow. However, from the present study, we find that the turbulent Prandtl number indeed depends on the molecular Prandtl number and this observation is based on the increasing value of ${\textit {Pr}}_t$ at the wall with respect to the scalars with ${\textit {Pr}}=4,6$, as shown in figure 13.

In order to verify the plausibility of the present observations at ${\textit {Pr}} = 4,6$, the obtained DNS data were compared with the DNS channel data reported by Alcántara-Ávila & Hoyas (Reference Alcántara-Ávila and Hoyas2021). Figure 15 shows the comparison of ${\textit {Pr}}_{t}$ at different molecular Prandtl numbers, where the present DNS data were at ${\textit {Re}}_{\theta } = 1070$, corresponding to ${\textit {Re}}_{\tau } = 395$, and the data from Alcántara-Ávila & Hoyas (Reference Alcántara-Ávila and Hoyas2021) were at ${\textit {Re}}_{\tau } = 500$. Despite these differences, the turbulent Prandtl numbers close to the wall are in good agreement, confirming that the ${\textit {Pr}}$-scaled wall-normal heat flux decreases with an increase in ${\textit {Pr}}$ for ${\textit {Pr}} \gtrapprox 4$, as stated by Alcántara-Ávila & Hoyas (Reference Alcántara-Ávila and Hoyas2021). Thus, the present observation confirms the constant behaviour of the turbulent Prandtl number very close to the wall and highlights its dependence on the molecular Prandtl number, which has been often ignored in turbulent heat-transfer calculations.

Figure 15. Comparison of turbulent Prandtl number against the channel DNS data at ${\textit {Re}}_\theta = 1070$ for different molecular Prandtl numbers. Present DNS for (——, blue) ${\textit {Pr}} = 1$, (——, orange) ${\textit {Pr}} = 2$, (——, green) ${\textit {Pr}} = 4$, (——, red) ${\textit {Pr}} = 6$, channel DNS data by Alcántara-Ávila & Hoyas (Reference Alcántara-Ávila and Hoyas2021) at ${\textit {Re}}_{\tau } = 497$ and ($\circ$, blue) ${\textit {Pr}} = 1$, ($\circ$, orange) ${\textit {Pr}} = 2$, ($\circ$, green) ${\textit {Pr}} = 4$, ($\circ$, red) ${\textit {Pr}} = 6$.

A brief discussion of the Reynolds-stress budget is provided in Appendix A.

3.5. Higher-order statistics, shear stress and heat flux

The higher-order statistics (specifically, the values of the third- and fourth-order moments of a quantity) provide information on the non-Gaussian behaviour of turbulence. The third- and fourth-order moments are also called the skewness and flatness, respectively, and for a statistically stationary variable $m$, they are defined as

(3.20)\begin{gather} \mathcal{S}(m) = \left.\left\langle m^3\right\rangle\right/\left\langle m^2\right\rangle^{3/2}, \end{gather}

(3.21)\begin{gather}\mathcal{F}(m) = \left.\left\langle m^4\right\rangle\right/\left\langle m^2\right\rangle^{2}. \end{gather}

The skewness and flatness of the streamwise and wall-normal velocity components are shown in figure 16 where, for the spanwise velocity, we observe a Gaussian behaviour in the overlap region, as expected. The higher-order moments obtained by Vreman & Kuerten (Reference Vreman and Kuerten2014) for the turbulent channel flow are compared with the present TBL results at ${\textit {Re}}_{\tau } = 180$ (corresponding to ${\textit {Re}}_\theta = 420$), showing a reasonable agreement. From figure 2, it should be highlighted that ${\textit {Re}}_\theta = 420$ is roughly at the beginning of the turbulent regime and the deviation of the profiles for $y^+>100$ as observed in figure 16 is due to the intermittent wake region which decays to the values for a Gaussian distribution as the free stream is approached. There are some differences observed in the overlap region where the skewness of $u$ is roughly $10\,\%$ higher than that reported by Vreman & Kuerten (Reference Vreman and Kuerten2014). Although the plots indicate an overall good agreement, some drastic differences are observed in the flatness of the wall-normal and spanwise velocity components close to the wall, where the turbulence is highly intermittent. The flatness of the wall-normal velocity component approaches a value of $\approx 29$ close to the wall for the data provided by Vreman & Kuerten (Reference Vreman and Kuerten2014) whereas, in the TBL, it converges to $\approx 20$. This is still lower than the value of $22$ obtained by Kim, Moin & Moser (Reference Kim, Moin and Moser1987). The skewness of the pressure fluctuations is roughly 20 % higher in the overlap region of the TBL compared with the channel. Further, the intermittency (flatness) of the pressure is higher than the velocity components as reported by Kim et al. (Reference Kim, Moin and Moser1987), whereas the data obtained with TBL show an offset of 15 % throughout and are lower than the flatness behaviour observed in the channel. The flatness factor for pressure fluctuation at the wall approaches a value of 4.5 in the present simulations as compared with the values of 4.7 (Li et al. Reference Li, Schlatter, Brandt and Henningson2009) and 4.9 (Schewe Reference Schewe1983), whereas a value of $\approx$5.2 is reported by Vreman & Kuerten (Reference Vreman and Kuerten2014).

Figure 16. Skewness and flatness for the streamwise, wall-normal and spanwise velocity components at ${\textit {Re}}_{\theta } = 420$ corresponding to ${\textit {Re}}_{\tau } = 180$. (——, blue) Present DNS, ($\circ$, blue) channel DNS data by Vreman & Kuerten (Reference Vreman and Kuerten2014) at ${\textit {Re}}_{\tau } = 180$.

The Reynolds shear stress is plotted in figure 17 and it is compared with the TBL data provided by Jiménez et al. (Reference Jiménez, Hoyas, Simens and Mizuno2010) and the channel DNS data by Kozuka et al. (Reference Kozuka, Seki and Kawamura2009). Even if there is a good agreement between the data in the inner region, there is a clear difference in the overlap region between the channel and the TBL. Such differences were also observed in the plots of the fluctuations in the streamwise and wall-normal directions shown in figure 8. The peak of the Reynolds shear stress is higher for the TBL compared with a channel flow, which indicates a higher momentum transfer by the fluctuating velocity field in the TBL. Looking at the energy budgets for shear stress component in both the channel and TBL corresponding to ${\textit {Re}}_\tau = 395$ (plots not shown here), we find that there is higher production in the TBL compared with channel flows (with uniform-heat-flux wall conditions) in the overlap region whereas the dissipation was of similar magnitude.

Figure 17. Inner-scaled Reynolds shear stress as a function of the wall-normal direction at ${\textit {Re}}_{\theta }=1070$. (——, blue) Present DNS, ($\circ$, blue) Jiménez et al. (Reference Jiménez, Hoyas, Simens and Mizuno2010) at ${\textit {Re}}_{\theta } = 1100$, ($\square$, blue) channel DNS data by Kozuka et al. (Reference Kozuka, Seki and Kawamura2009) at ${\textit {Re}}_{\tau } = 395$.

The turbulent streamwise heat flux shows a reasonable agreement with the channel data available at ${\textit {Pr}} = 1, 2$, as shown in figure 18, although the peaks for the channel-flow data are slightly higher and closer to the wall compared with the present observations. For ${\textit {Pr}} = 4$, however, the peak of the streamwise heat flux is higher in the channel data by Alcántara-Ávila et al. (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2018), since the heat flux is reported at a higher Reynolds number than our DNS. Furthermore, the comparison at the highest ${\textit {Pr}}$ in our simulation was at a slightly lower ${\textit {Pr}}$ than that in the channel and hence there is a difference in the slope of the streamwise heat flux close to the wall. Figure 19 shows the wall-normal heat flux for the different simulations. A clear difference can be observed in the outer region when comparing similar Prandtl numbers. On the other hand, the boundary-layer data and the channel data show a better agreement in the inner region.

Figure 18. Variation of streamwise heat flux along the wall-normal direction at ${\textit {Re}}_\theta = 1070$ for different passive scalars. Present DNS data for (——, blue) ${\textit {Pr}} = 1$, (——, orange) ${\textit {Pr}} = 2$, (——, green) ${\textit {Pr}} = 4$, (——, red) ${\textit {Pr}} = 6$, channel DNS data by Kozuka et al. (Reference Kozuka, Seki and Kawamura2009) at ${\textit {Re}}_{\tau } = 395$ and ($\circ$, blue) ${\textit {Pr}}=1$, ($\circ$, orange) ${\textit {Pr}}=2$, ($\circ$, red) ${\textit {Pr}}=7$, ($\circ$, green) channel DNS data by Alcántara-Ávila et al. (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2018) at ${\textit {Re}}_{\tau } = 500$ and ${\textit {Pr}} = 4$.

Figure 19. Variation of wall-normal heat flux along the wall-normal direction at ${\textit {Re}}_\theta = 1070$ for different passive scalars. Present DNS data for (——, blue) ${\textit {Pr}} = 1$, (——, orange) ${\textit {Pr}} = 2$, (——, green) ${\textit {Pr}} = 4$, (——, red) ${\textit {Pr}} = 6$, channel DNS data by Kozuka et al. (Reference Kozuka, Seki and Kawamura2009) at ${\textit {Re}}_{\tau } = 395$ and ($\circ$, blue) ${\textit {Pr}}=1$, ($\circ$, orange) ${\textit {Pr}}=2$, ($\circ$, red) ${\textit {Pr}}=7$, ($\circ$, green) channel DNS data by Alcántara-Ávila et al. (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2018) at ${\textit {Re}}_{\tau } = 500$ and ${\textit {Pr}} = 4$.

The correlation coefficients provide more information on the statistical association between the fields, and here the structure of the flow field and scalar fluctuations are analysed in terms of

(3.22)\begin{equation} \left. \begin{aligned} R_{u\theta_i} & = \frac{\left\langle u \theta_i\right\rangle}{u_{rms}\theta_{i\,{rms}}},\\ R_{v\theta_i} & = \frac{\left\langle v \theta_i\right\rangle}{v_{rms}\theta_{i\,{rms}}}. \end{aligned} \right\} \end{equation}

The correlation-coefficient plot in figure 20 corresponding to $u-\theta _i$ shows a strong correlation of streamwise velocity and scalar fluctuations at ${\textit {Pr}} = 1$ and it decreases with increase in Prandtl number. This is related to the similarity in the momentum and passive scalar transport by the turbulent eddies close to the wall. On the other hand, the $v-\theta _i$ correlation coefficient also shown in figure 20 exhibits an increasing trend with the Prandtl number. In the conductive sub-layer the correlation coefficients coincide for the various Prandtl numbers under study and then approach different values at the wall. This highlights a similar behaviour of the turbulent wall-normal momentum and passive scalar fields with a caveat, i.e. the differences are present very close to the wall for different Prandtl numbers.

Figure 20. Correlation coefficients: (a) $u-\theta _i$, (b) $v-\theta _i$. Present DNS data for scalars at (——, blue) ${\textit {Pr}} = 1$, (——, orange) ${\textit {Pr}} = 2$, (——, green) ${\textit {Pr}} = 4$, (——, red) ${\textit {Pr}} = 6$ at ${\textit {Re}}_\theta = 1070$.

3.6. Spanwise two-point correlations

The two-point correlations provide some quantitative information on the turbulent structures near the wall. For example, the streak spacing near the wall is of interest and can also be observed in an experimental setting. In order to identify the mean streak spacing, spanwise two-point correlations of the velocity components and passive scalars were obtained at five different positions along the streamwise direction at different wall-normal positions. Overall, the obtained results at ${\textit {Re}}_\theta = 830$ were compared with the data reported by Kim et al. (Reference Kim, Moin and Moser1987) and show good agreement (not depicted here). The obtained two-point correlations for different scalars at ${\textit {Re}}_\theta = 830$ are shown in figure 21 to assess the differences for varying ${\textit {Pr}}$. The two-point correlation becomes negative and reaches a minimum at an inner-scaled correlation length of $\delta z^+ \approx 50$. The length at which the minimum occurs provides an estimate of the half-mean separation between the streaks in the spanwise direction, i.e. $(\lambda _s^+/2)$. The streak spacing is plotted along the streamwise direction at a wall-normal position of $y^+ \approx 7$, as shown in figure 21. Overall, the velocity and scalar streak spacings increase with ${\textit {Re}}_\theta$, as reported in the works of Li et al. (Reference Li, Schlatter, Brandt and Henningson2009) and Schlatter et al. (Reference Schlatter, Örlü, Li, Brethouwer, Fransson, Johansson, Alfredsson and Henningson2009). Note that the correlations are available only at five streamwise locations (solid line for $U$ in figure 21 is for clarity), but the comparison with Li et al. (Reference Li, Schlatter, Brandt and Henningson2009) is in reasonably good agreement for the velocity and scalar streaks at ${\textit {Pr}} = 2$. The velocity streak spacing increases from 102 to 115 and appears to saturate for ${\textit {Re}}_\theta > 830$. Note that such a saturation of streak spacing was also pointed out by Li et al. (Reference Li, Schlatter, Brandt and Henningson2009) for ${\textit {Re}}_\theta > 1500$. From figure 21, we also observe that the streak spacing decreases with increasing ${\textit {Pr}}$ and that the streak spacings for the scalars at ${\textit {Pr}} = 4, 6$ are indistinguishable, although the rate of decay of the two-point correlations for the scalars is different. A higher grid resolution might be necessary to quantify the possible differences in the scalar streak spacing at higher ${\textit {Pr}}$.

Figure 21. Variation of mean spanwise streak spacing along the streamwise direction at $y^+ \approx 7$. (——, blue) $U$, ($\circ$, blue) $\varTheta _1$, ($\circ$, orange) $\varTheta _2$, ($\circ$, green) $\varTheta _3$, ($\circ$, red) $\varTheta _4$. Comparison data from Li et al. (Reference Li, Schlatter, Brandt and Henningson2009) for ($\times$, black) $U$, ($\square$, black) $\varTheta _2$. Inset: spanwise two-point correlation for (——, blue) $\varTheta _1$, (——, orange) $\varTheta _2$, (——, green) $\varTheta _3$, (——, red) $\varTheta _4$ at ${\textit {Re}}_\theta \approx 830$.

3.7. Scaling of wall-heat-flux fields

The wall-shear and heat-flux fields at different Prandtl numbers are normalized by subtracting the mean and dividing by the corresponding r.m.s. quantities.The normalization of a quantity $n$ is calculated as

(3.23)\begin{equation} \bar{n} = \frac{n-\left\langle n\right\rangle}{n_{rms}}. \end{equation}

For the results discussed, $n$ could indicate streamwise wall-shear stress, streamwise scalar flux, streamwise velocity or scalar concentrations at different Prandtl numbers. Figure 22 shows the instantaneous normalized streamwise wall-shear and heat-flux fields at different Prandtl numbers. The instantaneous normalized streamwise velocity and scalar fields are shown at $y^+=15$ (corresponding to $x/\delta _0^* = 500$) in figure 23. The normalized wall fields appear qualitatively similar and this result is confirmed by the distribution of the data in the streamwise shear and heat-flux fields obtained from 3700 samples, shown in figure 24. Although the distribution of data is different in the various fields, after normalization the distributions become identical. This indicates the uniformity in the distribution of the fluctuations of shear and heat flux when scaled with the corresponding r.m.s. quantities. This observation is useful for certain applications, for instance, in the prediction of fluctuating flow quantities from the wall, as discussed in the study by Guastoni et al. (Reference Guastoni, Balasubramanian, Güemes, Ianiro, Discetti, Schlatter, Azizpour and Vinuesa2022).

Figure 22. Instantaneous normalized streamwise wall-shear and heat-flux fields at different Prandtl numbers.

Figure 23. Instantaneous normalized streamwise velocity and passive scalar fields for different Prandtl numbers corresponding to $y^+ = 15$ at $x/\delta _0^* = 500$.

Figure 24. Probability density function of wall-shear and heat-flux fields at different Prandtl numbers. (a) Distribution of absolute data in the fields with a bin size of 1 unit, (b) distribution of the normalized data in the fields with a bin size of 0.1 unit, (——, blue) ${\partial U}/{\partial y}$ or $\overline {{\partial U}/{\partial y}}$, (——, orange) ${\partial \varTheta _{1}}/{\partial y}$ or $\overline {{\partial \varTheta _{1}}/{\partial y}}$, (——, green) ${{\partial \varTheta _{2}}/{\partial y}}$ or $\overline {{\partial \varTheta _{2}}/{\partial y}}$, (——, red) ${{\partial \varTheta _{3}}/{\partial y}}$ or $\overline {{\partial \varTheta _{3}}/{\partial y}}$, (——, purple) ${{\partial \varTheta _{4}}/{\partial y}}$ or $\overline {{\partial \varTheta _{4}}/{\partial y}}$.