2.1. Flow conditions

Table 1 outlines the free-stream and wall conditions of the novel DNS case, with $M_{\infty } = 4.9$ and $T_w/T_r = 0.60$ . With recent increases in available computational hardware, the Reynolds number range for the present DNS is $1000 \lesssim {Re}_\tau \lesssim 2500$ , distributed over continuous streamwise domains. The present case maintains the Mach number ( $M_\infty = 4.9$ ) of a case that was reported by Huang et al. (Reference Huang, Duan and Choudhari2022) as a numerical replication of the high-speed wind tunnel experiment performed at the National Aerothermochemistry Laboratory (NAL) at Texas A&M University. However, unlike the experimental campaign that featured a quasi-adiabatic wall, the current DNS has a colder wall temperature of $T_w/T_r = 0.60$ . To gauge the heat-transfer effects in compressible TBLs, Gibis et al. (Reference Gibis, Sciacovelli, Kloker and Wenzel2024) developed a novel classification scheme based on the Eckert number $Ec$ (being equivalent to diabatic parameter $\Theta$ ), which effectively categorises compressible cold-wall TBLs into weakly/moderately/strongly cooled cases. For the given $M_\infty$ and $T_w/T_r$ , the Eckert number and diabatic parameter are $-1/Ec\approx -0.223$ and $\Theta \approx 0.5$ , respectively. Figure 1 shows that the current DNS falls into the moderately cooled regime over the specified Reynolds number range. For comparison, the DNS of a $M_\infty =5.86$ TBL with ${Re}_\tau \simeq 2000$ by Cogo et al. (Reference Cogo, Salvadore, Picano and Bernardini2022) lies in the weakly cooled regime. Therefore, the current DNS complements existing studies by uniquely allowing for insight into the effect of Reynolds number in a different wall-cooling regime than before.

Table 1. Free-stream and wall-temperature conditions for the novel M5Tw060 DNS case: $T_r$ is the recovery temperature, $T_r = T_\infty [ 1 + r (\gamma - 1) M_\infty ^2/2 ]$ with $r=0.89$ ; $Ec = (\gamma - 1) M_\infty ^2 T_\infty /(T_r-T_w)$ is the Eckert number; $\Theta = (T_w - T_\infty )/(T_r - T_\infty )$ is the diabatic parameter. The various Reynolds numbers are defined as: ${Re}_u = \rho _\infty U_\infty /{\rm \mu} _\infty$ ; ${Re}_\theta = \rho _\infty U_\infty \theta /{\rm \mu} _\infty$ ; ${\rm {Re}}_\tau = \rho _w u_\tau \delta /{\rm \mu} _w$ ; ${Re}^*_\tau = \rho _\delta \sqrt {\tau _w/\rho _\delta } \delta /{\rm \mu} _\delta$ ; ${Re}_{\delta 2} = \rho _\infty U_\infty \theta /{\rm \mu} _w$ . The subscripts $\infty , w, \delta$ denote values in the free stream, at the wall and at the boundary-layer edge ( $z=\delta$ ), respectively. The specified range of Reynolds numbers corresponds to the `useful’ portion of the computational domain (i.e. from downstream of the inflow adjustment zone to the end of the computational domain).

Figure 1. Regime diagram (adapted from Gibis et al. Reference Gibis, Sciacovelli, Kloker and Wenzel2024) for various cooled wall cases, including the current study (indicated by green circles). Data are plotted versus $-1/E_c$ (left axis) and equivalently as the diabatic parameter $\Theta$ (right axis).

2.2. Governing equations and numerical methods

The DNS code used in this work solves the three-dimensional compressible Navier–Stokes equations in conservative form. The working fluid is air, modelled as a perfect gas subject to the usual constitutive relations for a Newtonian fluid. To compute the temperature--viscosity relationship, Sutherland’s law is used and the thermal conductivity from $k = {\rm {\rm \mu}} c_p/Pr$ , with $c_p$ as the heat capacity at constant pressure and $Pr = 0.71$ . A seventh-order weighted essentially non-oscillatory (WENO) scheme is used to discretise the inviscid fluxes. To reduce numerical dissipation, the current scheme is optimised using limiters (Taylor, Wu & Martín Reference Taylor, Wu and Martín2007; Wu & Martin Reference Wu and Martin2007), compared with the original WENO introduced by Jiang & Shu (Reference Jiang and Shu1996). The viscous fluxes are discretised using a fourth-order central difference scheme, and the time integration is performed with a third-order, low-storage Runge–Kutta scheme (Williamson Reference Williamson1980). Fundamentals of the numerical method are described in Martín (Reference Martín2007). The DNS code has been validated for supersonic and hypersonic TBLs (Wu & Martin Reference Wu and Martin2007; Duan et al. Reference Duan, Beekman and Martin2010, Reference Duan, Beekman and Martin2011; Duan & Martin Reference Duan and Martin2011; Duan et al. Reference Duan, Choudhari and Wu2014, Reference Duan, Choudhari and Zhang2016; Zhang et al. Reference Zhang, Duan and Choudhari2017, Reference Zhang, Duan and Choudhari2018; Huang et al. Reference Huang, Duan and Choudhari2022).

Figure 2. Computational domain and simulation set-up for case M5Tw060. The instantaneous flow is shown using isosurfaces of the magnitude of the density gradient, $|\nabla \rho |\delta _i/\rho _\infty \approx 0.98$ , coloured by the streamwise velocity component (with levels from 0 to $U_\infty$ , blue to red). The inflow boundary-layer thickness of box 1 DNS is $\delta _i = 2.65$ mm. The planes labelled as BL1, BL2 and BL3 are the analysis planes outlined in table 3.

2.3. Computational domain and simulation set-up

Figure 2 shows the computational domains used for the novel DNS case, referred throughout as M5Tw060. A series of streamwise overlapping boxes are used to cover a total distance of $L_x = 400.5\delta _i$ , $L_y = 9.96\delta _i$ , $L_z = 53.2\delta _i$ in the streamwise ( $x$ ), spanwise ( $y$ ) and wall-normal ( $z$ ) directions, respectively. The incoming boundary-layer thickness to the first box is $\delta _i = 2.65$ mm. A modified recycling--rescaling method (Duan, Choudhari & Wu Reference Duan, Choudhari and Wu2014) is used to bypass transition and generate turbulent structures in the first box. The recycling plane is placed at $26.4\delta _i$ downstream of the inflow plane (i.e. $(x_{rec}-x_i)/\delta _i=26.4$ ). All statistical analyses in this work are performed sufficiently downstream of the inflow plane ( $\gt 53\delta _i$ ) to avoid any unrealistic effects caused by the recycling--rescaling method. The subsequent downstream boxes serve to provide streamwise distance for the TBL to evolve. To mitigate any bias due to each overlapping box, spanwise/wall-normal planes are sampled adequately upstream of the end of each box and passed directly onto the next box. The resulting TBL is continuous across the boxes due to the streamwise overlap. At the upper and outflow boundaries, unsteady non-reflecting boundary conditions based on Thompson (Reference Thompson1987) are used, and periodic boundary conditions are applied in the spanwise direction. At the wall, no-slip conditions are imposed for the three velocity components, and an isothermal condition is used for the temperature with $T_w/T_r \simeq 0.6$ . A similar numerical set-up and boundary conditions were used in the previous DNS study by Huang et al. (Reference Huang, Duan and Choudhari2022).

For all DNS boxes, constant grid spacing is maintained in the streamwise and spanwise directions, with non-dimensional grid spacings of $\Delta x^+ = 6.5$ and $\Delta y^+ = 5.1$ , respectively. These spacings are non-dimensionalised using the boundary-layer thickness at the furthest downstream location (BL3) where statistics are sampled, as outlined in table 3. The wall-normal grid distributions are clustered near the wall with $\Delta z^+ = 0.52$ and are smoothly stretched to a larger spacing value of $\Delta z^+ = 8.3$ near the boundary-layer edge. An identical time step of $\Delta t^+ = 2.75 \times 10^{-2}$ is used for each DNS box. Both the instantaneous flow fields and the averaged turbulence statistics in the overlapping regions between each of the two neighbouring boxes are carefully examined to ensure a seamless connection of the computed flow fields across all boxes. The same multiple-domain method has been successfully employed and validated in our previous studies (Zhang et al. Reference Zhang, Duan and Choudhari2018; Huang et al. Reference Huang, Duan and Choudhari2022; Nicholson, Duan & Bisek Reference Nicholson, Duan and Bisek2024).

Table 2. Summary of parameters for the Mach $4.9$ DNS dataset. Here $L_x$ , $L_y$ and $L_z$ are the domain sizes in the streamwise, spanwise and wall-normal directions, respectively; $N_x$ , $N_y$ and $N_z$ are the grid dimensions; $\Delta x^+$ and $\Delta y^+$ are the uniform grid spacing in the streamwise and spanwise directions; $\Delta z^+$ denotes the wall-normal spacing at the first grid away from the wall and near the boundary-layer edge; $L_{to} = U_\infty \delta _i / u_{\tau ,i}$ is one turnover length of the largest eddy at inflow, where $\delta _i = 2.65$ mm is the inflow boundary-layer thickness and $u_{\tau ,i} = 36.9$ m s⁻¹ is the inflow friction velocity. All grid spacings are normalised by the viscous scale at the farthest downstream station selected for statistical analysis as listed in table 3.

Table 3. Boundary-layer properties at the DNS station $x_a$ selected for analysis. Here $x_i=0.151654$ m denotes the streamwise coordinate at the inflow plane. The boundary-layer thickness $\delta$ is defined as the wall-normal distance from the wall to the location where $\overline {u}=0.99U_{\infty }$ ; $H_{12}= \delta ^*/\theta$ is the shape factor; $u_\tau =\sqrt {\tau _w/\rho _w}$ is the friction velocity; $z_\tau =\nu _w/u_\tau$ is the viscous length. Reynolds numbers are defined in Table 1.

The turbulent statistics presented in this paper are first averaged in the homogeneous spanwise direction, followed by a temporal average over $549$ snapshots spanning a non-dimensional time interval of $T_f u_{\tau , i}/ \delta _i = 38.2$ , where $T_f$ is the total time considered for collecting the flow statistics. To improve statistical convergence, particularly for the streamwise evolution and the profile plots presented in § 3, the statistics are further averaged within a streamwise window of $[x_a - 0.5\delta , x_a + 0.5\delta ]$ , where $x_a$ is the reference streamwise location and $\delta$ is the local boundary-layer thickness. This approach is widely used to ensure statistical convergence (Pirozzoli Reference Pirozzoli2011; Duan, Choudhari & Zhang Reference Duan, Choudhari and Zhang2016; Huang et al. Reference Huang, Duan and Choudhari2022). Statistical convergence is validated by varying the number of snapshots significantly and ensuring that the differences in the computed statistics across various intervals remain minimal. In this paper the Reynolds average of a state variable $f$ is denoted as $\bar {f}$ and the Favre average of that is represented by $\widetilde {f}$ , where $\widetilde {f}=\overline {\rho f}/\bar {\rho }$ . The fluctuations around these averages are denoted by $f^\prime$ and $f^{\prime \prime }$ , respectively.

The power spectra presented in §§ 4 and 5 are computed using the Welch method (Welch Reference Welch1967), where the data are divided into four segments, each spanning $411.2\delta _i/U_\infty$ with a $50\,\%$ overlap between consecutive segments. A standard Hanning window is used to minimise spectral leakage prior to performing the fast Fourier transform. The sampling frequency is approximately $67U_\infty /\delta _i$ (corresponding to $\Delta t_{sample}^+ = 2.75 \times 10^{-1}$ ), while the total signal length is around $1028\delta _i/U_\infty$ .

Additional information about the validation of the DNS, including an assessment of the domain size as well as comparisons with available data in the literature, is given in the Appendix.

3.1. Skin friction and surface heat flux

First, we use the DNS data to examine the effect of Reynolds number on skin friction and surface heat flux of the TBL.

The van Driest II transformation (Van Driest Reference Van Driest1956) and the Spalding & Chi (Reference Spalding and Chi1964) transformations are commonly used to transform the skin friction associated with a compressible boundary layer into that of an equivalent, incompressible boundary layer. Figure 4 compares the performance of van Driest II and Spalding and Chi transformations for the current DNS case (case M5Tw060). Additional cases from Huang et al. (Reference Huang, Duan and Choudhari2022) are included, which feature various Mach numbers and wall temperatures ranging from $2.5 \leqslant M_\infty \leqslant 14$ and $0.20 \leqslant T_w/T_r \leqslant 1.0$ . After the skin friction associated with each high-speed TBL is transformed into that of an equivalent, incompressible boundary layer, it is compared against the incompressible correlations of Karman–Schoenherr (Roy & Blottner Reference Roy and Blottner2006), Smits (Smits, Matheson & Joubert Reference Smits, Matheson and Joubert1983) and Coles–Fernholtz (Nagib, Chauhan & Monkewitz Reference Nagib, Chauhan and Monkewitz2007). The detailed formulations for the van Driest II and Spalding and Chi transformations, along with those for the aforementioned incompressible friction correlations, are listed in Huang et al. (Reference Huang, Duan and Choudhari2022). The current DNS corroborates with the previous DNS cases and features almost equally good agreement with the incompressible correlations of Smits and Coles–Fernholtz after the van Driest II transformation is performed, with no indication that this trend will not continue at even higher Reynolds numbers. The Spalding and Chi transformation tends to slightly underpredict $C_f$ by approximately $10\,\%\!-\!20\,\%$ at lower Reynolds numbers. Interestingly, for the case M5Tw060, the discrepancy between the Spalding and Chi transformed skin friction values of the DNS and the incompressible correlations decreases gradually with increasing Reynolds number. Such a finding is consistent with the reports of Goyne, Stalker & Paull (Reference Goyne, Stalker and Paull2003) and Huang et al. (Reference Huang, Duan and Choudhari2022).

Figure 4. Transformed skin friction coefficient ( $C_{f,i}=F_c C_f$ ) versus Reynolds numbers ( ${Re}_{\theta ,i}=F_\theta {Re}_{\theta }$ ) based on (a) the van Driest (VD) II theory and (b) the Spalding and Chi (SC) theory, wherein the black solid line, dash-dotted line and dashed line denote the incompressible correlations of Karman–Schoenherr (Roy & Blottner Reference Roy and Blottner2006), Coles–Fernholtz (Nagib et al. Reference Nagib, Chauhan and Monkewitz2007) and Smits (Smits et al. Reference Smits, Matheson and Joubert1983), respectively.

Figure 5. Comparison of (a) skin friction coefficient $C_f$ and (b) Stanton number $C_h$ versus Reynolds number in comparison with the van Driest (VD) II theory and the predictive models from Kumar & Larsson (Reference Kumar and Larsson2022) and Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2024).

Kumar & Larsson (Reference Kumar and Larsson2022) and Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2024) recently proposed models for predicting drag and heat transfer in compressible high-speed flows, utilising a compressible velocity scaling law that inverse transforms the velocity profile of incompressible flows and the temperature--velocity relationship of Zhang et al. (Reference Zhang, Bi, Hussain and She2014). Figures 5(a) and 5(b) compare DNS and model predictions for the compressible skin friction coefficient $C_f$ and Stanton number $C_h$ as a function of the Reynolds numbers, respectively. Here, $C_f=2\tau _w/\rho _\infty U_\infty ^2$ is the skin friction coefficient and $C_h=q_w/[\rho _\infty U_\infty c_p (T_r-T_w)]$ is the Stanton number. The models of Kumar & Larsson (Reference Kumar and Larsson2022) and Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2024) demonstrate commendable accuracy in predicting $C_f$ and $C_h$ within $\pm 4\,\%$ , particularly over the higher-Reynolds-number range ( ${Re}_\theta \gtrsim 18\,000$ or ${Re}_\tau \gtrsim 1500$ ) for this case. Such an excellent prediction confirms the reliability of the models under high-Reynolds-number conditions.

The Reynolds analogy factor $R_{af}=2C_h/C_f$ is often used to obtain the surface heat flux from the skin friction coefficient (Roy & Blottner Reference Roy and Blottner2006). Typically, the Reynolds analogy factor is assumed to be $R_{af} = Pr^{-2/3}$ ( ${\approx} 1.256$ for $Pr = 0.71$ ), or, from the reports of Chi & Spalding (Reference Chi and Spalding1966) and Hopkins & Inouye (Reference Hopkins and Inouye1971), $R_{af} = 1.16$ . The recent DNS study of hypersonic TBLs with cooled walls by Huang et al. (Reference Huang, Duan and Choudhari2022) reported values of the Reynolds analogy factor to be within the narrow range of $1.1 \lt R_{af} \lt 1.2$ up to ${Re}_\tau \approx 1200$ . Figure 6(a) shows the evolution of $R_{af}$ versus friction Reynolds number ${Re}_\tau$ for the present DNS case, along with some of the DNS cases in Huang et al. (Reference Huang, Duan and Choudhari2022). For the present DNS, the Reynolds analogy factor follows a relatively constant trend at $R_{af} \approx 1.12$ , which is consistent with the previous cases (i.e. $1.1 \lt R_{af} \lt 1.2$ ), whereas the lower-Reynolds-number cases with highly cooled walls depict a slightly negative slope, tending towards the constant value of $R_{af} \approx 1.16$ up to ${Re}_\tau \approx 1200$ . This different behaviour of $R_{af}$ for the present DNS with respect to previous cases with highly cooled walls may be attributed to differences in Reynolds number ${Re}_\tau$ and wall-cooling rates, considering that $R_{af}$ assumes smaller values at larger $\Theta$ as shown by the data from Cogo et al. (Reference Cogo, Baù, Chinappi, Bernardini and Picano2023), and the variation in the Reynolds number can also have an impact on the heat-transfer effects at a fixed $\Theta$ or $-1/Ec$ (Gibis et al. Reference Gibis, Sciacovelli, Kloker and Wenzel2024). Figure 6(b) further plots the computed Reynolds analogy factor ( $R_{af} Pr$ ) as a function of the diabatic parameter $\Theta$ in comparison with those reported in Cogo et al. (Reference Cogo, Baù, Chinappi, Bernardini and Picano2023). The predicted Reynolds analogy factor by the current DNS is consistent with those reported in Cogo et al. (Reference Cogo, Baù, Chinappi, Bernardini and Picano2023) with a similar diabatic parameter of $\Theta = 0.5$ . It also shows good agreement with the empirical correlation proposed by Zhang et al. (Reference Zhang, Bi, Hussain and She2014), which bounds $R_{af} Pr$ within the range of $0.8 \pm 0.03$ .

Figure 6. (a) The Reynolds analogy factor $R_{af}=2C_h/C_f$ as a function of friction Reynolds number ${Re}_\tau$ where the horizontal dashed and dashed dot lines denote constant values of $1.16$ and $1.12$ , respectively. (b) Variation of $(2C_h/C_f) Pr$ with the diabatic parameter, $\Theta = (T_w - T_\infty )/(T_r - T_\infty )$ . DNS results from Cogo et al. (Reference Cogo, Baù, Chinappi, Bernardini and Picano2023) are included for comparison. The shaded region represents the range of $0.8 \pm 0.03$ according to Zhang et al. (Reference Zhang, Bi, Hussain and She2014).

3.2. Mean velocity and temperature

Next, the mean velocity profiles as a function of Reynolds number are investigated, focusing on assessing compressibility transformations that map the mean velocity profile in a compressible turbulent boundary to an equivalent incompressible profile. At the centre of the compressibility transformations are Morkovin’s hypothesis (Morkovin Reference Morkovin1962) and the van Driest transformation (Van Driest Reference Van Driest1951). While the van Driest transformation is a classical velocity scaling law that is known to successfully collapse the inner and log layer of incompressible TBLs to the law of the wall for adiabatic boundary layers, it is unable to recover the incompressible law of the wall in compressible flows with non-adiabatic wall conditions (i.e. flows with surface heat transfer) (Duan et al. Reference Duan, Beekman and Martin2010; Zhang et al. Reference Zhang, Duan and Choudhari2018). More recent transformations that have shown to better account for the surface heat-transfer effect include those of Trettel & Larsson (Reference Trettel and Larsson2016), Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020) and Griffin, Fu & Moin (Reference Griffin, Fu and Moin2021). Cogo et al. (Reference Cogo, Salvadore, Picano and Bernardini2022) recently applied the aforementioned transformations to their novel dataset up to ${Re}_\tau \approx 2000$ at $M_\infty = 5.86$ with a weakly cooled wall, and concluded that the more recent transformations such as those by Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020) and Griffin et al. (Reference Griffin, Fu and Moin2021) appear to be more efficient in capturing the effect of compressibility and collapsing the velocity profiles. Similar conclusions are drawn in the present work with a higher Reynolds number of ${Re}_\tau \approx 2500$ and a moderately cooled wall. The detailed comparison of the current DNS with the mean velocity transformations of Trettel & Larsson (Reference Trettel and Larsson2016), Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020) and Griffin et al. (Reference Griffin, Fu and Moin2021) are outlined in the Appendix (figure 35).

Recently, Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2023) introduced a novel velocity transformation (hereafter referred to as the Hasan-Larsson-Pirozzoli-Pecnik (HLPP) transformation) that accounts for intrinsic compressibility effects by identifying the friction Mach number $M_\tau$ as the key parameter and adjusting the near-wall damping function accordingly. The mean velocity profile is given by

(3.1a) \begin{align} u_{{{\rm HLPP}}}^{+} &= \int _{0}^{{u}^{+}} \left (\frac {1 + {\rm {\rm \mu}} ^c_t/\overline {{\rm \mu} }}{1 + {\rm \mu} ^i_t/{\rm \mu} _w}\right ) \left (1 - \frac {z}{z_{\tau }^{*}}\frac {{\rm d}z_\tau ^{*}}{{\rm d}z}\right ) \sqrt {\frac {\overline {\rho }}{\rho _{w}}} {\rm d}{u}^{+} \end{align}

(3.1b) \begin{align} & = \int _{0}^{{u}^{+}} \left (\frac {1 + \kappa z^{*} D^{c}}{1 + \kappa z^{*} D^{i}}\right ) \left (1 - \frac {z}{z_{\tau }^{*}}\frac {{\rm d}z_\tau ^{*}}{{\rm d}z}\right ) \sqrt {\frac {\overline {\rho }}{\rho _{w}}} {\rm d}{u}^{+}, \end{align}

where $\kappa =0.41$ is the von Kármán constant, and $D^i$ and $D^c$ represent the damping function for incompressible and compressible flows, respectively. These functions are defined as

(3.2a) \begin{align} D^{i} &= \left [1 - \exp \left (\frac {-z^{*}}{A^{+}}\right )\right ]^{2}, \end{align}

(3.2b) \begin{align} D^{c} &= \left [1 - \exp \left (\frac {-z^{*}}{A^{+} + f(M_{\tau })}\right )\right ]^{2}, \end{align}

with $A^+=17$ and $f(M_\tau )=19.3M_\tau$ .

Figure 7(a) presents the mean velocity profiles transformed using the recently proposed HLPP transformation by Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2023). The HLPP transformation effectively captures the effects of compressibility and achieves a very good collapse of velocity profiles across the Reynolds number range, giving the von Kármán constant $\kappa = 0.41$ and a corresponding log-law intercept of $5.2$ that is consistent with the classical incompressible log-law intercept. It provides a more accurate representation of the logarithmic region by effectively eliminating the log-law shift that was observed in the Trettel and Larsson transformation (Zhang et al. Reference Zhang, Duan and Choudhari2018).

Figure 7. (a) Mean velocity profiles transformed using the method proposed by Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2023). Incompressible DNS data from Sillero, Jiménez & Moser (Reference Sillero, Jiménez and Moser2013) at a comparable friction Reynolds number is also plotted for comparison. (b) The diagnostic function $I_c = \hat {Z}({\rm d} u_c^+/{\rm d} \hat {Z})$ of HLPP transformed velocity profiles with $u_c^+ = u_{\textit{HLPP}}^+$ and $\hat {Z} = z^*$ .

To quantitatively investigate the extent of the log region, a log-law diagnostic function is used, similar to that of Pirozzoli & Bernardini (Reference Pirozzoli and Bernardini2013) and Huang et al. (Reference Huang, Duan and Choudhari2022):

(3.3) \begin{equation} I_c = \hat {Z}\bigg(\frac{{\rm d} u_c^+}{{\rm d} \hat {Z}}\bigg). \end{equation}

Here $u_c^+$ and $\hat {Z}$ are the transformed mean velocity and non-dimensional wall distance, respectively, for a given compressibility transformation. In this definition, a constant value of $I_c$ indicates perfectly logarithmic behaviour. Figure 7(b) shows that the recent HLPP transformation exhibits a nearly constant $I_c$ over an wide range ( $z^* \approx 50{-}350$ ), and the region of constancy is considerably extended at the highest Reynolds number. The extended region of constancy for the diagnostic function is also seen for several other mean velocity transformations such as those of Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020) and Griffin et al. (Reference Griffin, Fu and Moin2021). For completeness, the diagnostic functions for these other transformations are outlined in the Appendix (figure 36).

One of the most dominant consequence of wall cooling is the near-wall sign change in the wall-normal temperature gradient compared with adiabatic conditions, which causes the near-wall regions with $\partial \widetilde {T}/\partial {z}\gt 0$ and negative turbulent heat flux $\overline {\rho T^{\prime\prime} w''}\lt 0$ (Duan et al. Reference Duan, Beekman and Martin2010; Gibis et al. Reference Gibis, Sciacovelli, Kloker and Wenzel2024). In a moderately cooled TBL in particular, the near-wall temperature peak, located at $\partial \widetilde {T}/\partial {z}=0$ , lies within the buffer layer ( $5 \lesssim z^* \lesssim 30$ ) where there is strong turbulence production, while for a weakly cooled TBL, $z^*_{\partial \widetilde {T}/\partial {z}=0}$ remains in the viscous sublayer ( $z^* \lesssim 5$ ) where the behaviour of the boundary layer is dominated by viscosity rather than turbulence. Figure 8 shows that the temperature peak or $z^*_{\partial \widetilde {T}/\partial {z}=0}$ is located within $5 \lesssim z^* \lesssim 30$ for the current DNS case, confirming that the DNS lies in the moderately cooled wall regime. The location of $z^*_{\partial \widetilde {T}/\partial {z}=0}$ is shown to shift to larger $z^*$ values as ${Re}_\tau$ increases, indicating stronger wall cooling at higher Reynolds numbers. Such a trend is consistent with the regime diagram of Gibis et al. (Reference Gibis, Sciacovelli, Kloker and Wenzel2024). Additionally, the theoretical model of Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2024) performed reasonably well in predicting the mean temperature profile as well as the outward shift of $z^*_{\partial \widetilde {T}/\partial {z}=0}$ with ${Re}_\tau$ .

Figure 8. Profiles of the (a) mean temperature $\widetilde {T}/T_w$ , (b) premultiplied temperature gradient $z^* (\partial \widetilde {T}^+/ \partial {z^*}) (\bar {\rho }/{\rho }_w)$ and (c) the peak temperature location $z^*_{\partial \widetilde {T}/\partial {z}=0}$ as a function of Reynolds number. In (a) and (b) the shaded region indicates the region of $5\leqslant z^*\leqslant 30$ , and the DNS data are compared with the theoretical prediction of Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2024). In (b) the temperature gradient is presented in wall units where $\widetilde {T}^+=\widetilde {T}/\left(\gamma T_w u_\tau /\sqrt {(\gamma R {\widetilde {T}}^{^{^{^{^{\!\!\!}}}}})}\right)$ . In (c) the previous DNS of Huang et al. (Reference Huang, Duan and Choudhari2022), Cogo et al. (Reference Cogo, Baù, Chinappi, Bernardini and Picano2023) and Gibis et al. (Reference Gibis, Sciacovelli, Kloker and Wenzel2024) are included for comparison.

Figure 9. Profiles of Reynolds stress components $\overline {\rho u^{\prime\prime}_i u^{\prime\prime}_j}/\tau _w$ as a function of Reynolds number. In (a), black dashed lines denote $\overline {\rho u''u''}/\tau _w=2.39-1.03\log (z/\delta )$ (Pirozzoli & Bernardini Reference Pirozzoli and Bernardini2011). In (b) the lines represent $\overline {\rho v''v''}/\tau _w=B_2-0.27\log (z/\delta )$ with $B_2=1.3$ (dashed line) and $B_2=1.5$ (dash-dotted line) (Baidya et al. Reference Baidya, Philip, Hutchins, Monty and Marusic2021).

3.3. Reynolds stresses

Figure 9 shows the normal and shear components of the Reynolds stress tensor $\overline {\rho u^{\prime\prime}_i u^{\prime\prime}_j}/\tau _w$ as a function of Reynolds number, where Morkovin’s hypothesis and semi-local scaling are applied to account for the compressibility effect. In all of the components of the Reynolds stress tensor investigated herein, the effect of Reynolds number is apparent via heightened turbulent fluctuations and intensities with increasing ${Re}_\tau$ , which are predominantly attributed to large-scale turbulent motions that convey their energetic footprint onto the smaller, near-wall scales. Based on the attached-eddy hypothesis (Townsend Reference Townsend1976), Perry & Li (Reference Perry and Li1990) demonstrated through experimentation that the variance of the velocity fluctuations should follow a logarithmic scaling as (Pirozzoli & Bernardini Reference Pirozzoli and Bernardini2011)

(3.4) \begin{equation} \frac {\overline {u_i^{{\prime}2}}}{u_\tau ^2}=B_i - A_i \log \left(\frac{z}{\delta} \right) -V(z^+), \end{equation}

where $B_i$ and $A_i$ are constants with $i = 1$ and $i =2$ denoting the streamwise and spanwise directions, respectively, and $V(z^+)$ accounts for viscous corrections and is set to zero. In the streamwise direction the constants are determined to be $B_1=2.39$ , $A_1=1.03$ (Pirozzoli & Bernardini Reference Pirozzoli and Bernardini2011). For the spanwise direction, a recent study of incompressible boundary layers by Baidya et al. (Reference Baidya, Philip, Hutchins, Monty and Marusic2021) proposed a constant value of $A_2=0.27$ and an increasing value of $B_2$ with ${Re}_\tau$ before $B_2$ asymptotes to approximately $1.5$ at ${Re}_\tau \approx 10\,000$ . Figures 9(a) and 9(b) demonstrate the apparent formation of logarithmic scaling with ${Re}_\tau$ for the streamwise and spanwise Reynolds stresses, respectively. A better approximation of the logarithmic law is seen for the spanwise velocity component ( $\overline {\rho v''v''}$ ) than for the streamwise component ( $\overline {\rho u''u''}$ ), which is consistent with the Mach 2 data by Pirozzoli & Bernardini (Reference Pirozzoli and Bernardini2011). Also, consistent with the incompressible data reported in Baidya et al. (Reference Baidya, Philip, Hutchins, Monty and Marusic2021), the value of $B_2$ increases with ${Re}_\tau$ . At ${Re}_\tau = 2480$ , $\overline {\rho v''v''}$ for the current data closely displays logarithmic scaling with $B_2=1.5$ , which is similar to the value determined for incompressible boundary layers at ${Re}_\tau \approx 10\,000$ (Baidya et al. Reference Baidya, Philip, Hutchins, Monty and Marusic2021). Such a trend suggests that the asymptotic high-Reynolds-number limit to a constant $B_2$ is close to be reached for the current flow, and the Morkovin’s hypothesis and semi-local scaling have successfully accounted for the effect of compressibility on the spanwise component of Reynolds stress. The presence of a noticeable bump in the outer region of $\overline {\rho v''v''}$ is also clear with increasing Reynolds number, which is likely a footprint of the large-scale motions in the wake region of the boundary layer (Pirozzoli & Bernardini Reference Pirozzoli and Bernardini2011). With regard to the wall-normal and shear components of the Reynolds stress (figures 9 c and 9 d, respectively), increasing the Reynolds number gives rise to a plateau of peak magnitude that extends for approximately a decade of wall-normal values at the highest Reynolds number, and such a plateau corresponds to the formation of an equilibrium layer in the boundary layer (Pirozzoli & Bernardini Reference Pirozzoli and Bernardini2011, Reference Pirozzoli and Bernardini2013; Cogo et al. Reference Cogo, Salvadore, Picano and Bernardini2022). Similar to the spanwise Reynolds stress, the large-scale dynamics in the outer most region of the boundary layer leaves a noticeable bump in the distribution of wall-normal Reynolds stress in the outer layer.

Figure 10. (a) The ratio of eddy viscosity to molecular viscosity ${\rm \mu} _t / \overline {{\rm \mu} }$ and (b) the total shear stress $[{\overline {{\rm \mu} } ({\rm d}\widetilde {u}/{\rm d}z)}-\overline {\rho u''w''}]/\tau _w$ as a function of Reynolds number. In (a) the square and diamond symbols denote the theoretical models of Johnson & King (Reference Johnson and King1985) and Huang et al. (Reference Huang, Coleman, Spalart and Yang2023), respectively.

Figure 10(a) presents the eddy viscosity ${\rm \mu} _t$ as a function of Reynolds number. Here, the eddy viscosity is computed from the DNS results using the following expression:

(3.5) \begin{equation} {\rm \mu} _t = -\frac {\overline {\rho u'' w''}}{{\rm d}\widetilde {u}/{\rm d} z}. \end{equation}

Here $\overline {\rho u'' w''}$ is the Reynolds shear stress and ${\rm d} \widetilde {u}/{\rm d} z$ denotes the primary mean strain rate. For comparison, we also include ${\rm \mu} _t$ predicted by the classic model of Johnson & King (Reference Johnson and King1985) and the recent model of Huang et al. (Reference Huang, Coleman, Spalart and Yang2023), which are defined as

(3.6a) \begin{align} \text{Johnson-King: } \quad {\rm \mu} _t &= \kappa z^+ D \overline {{\rm \mu} }, \quad \quad \text{with } D = [1 - \exp (-z^+ / A^+)]^2, \end{align}

(3.6b) \begin{align} \text{Huang } et\, al.: \quad \quad {\rm \mu} _t &= \kappa z^* D^* \overline {{\rm \mu} }, \quad \quad \text{with } D^* = [ 1 - \exp ( -{z^*}/{A^*})]^2. \end{align}

Here, $\kappa = 0.41$ is the von Kármán constant, $\overline {{\rm \mu} }$ is the mean molecular viscosity calculated based on Sutherland’s law, $D$ and $D^*$ are damping functions based on the inner ( $z^+$ ) and semi-local ( $z^*$ ) scales, respectively, with model constants of $A^+ = 17$ and $A^*=17$ . As expected, the theoretical model of Huang et al. (Reference Huang, Coleman, Spalart and Yang2023), which is defined using a semi-local wall distance $z^*$ , provides a much improved prediction of ${\rm \mu} _t / \overline {{\rm \mu} }$ in the inner region of the boundary layer than the classic model of Johnson & King (Reference Johnson and King1985). The wall-normal extent with good model performance increases with ${Re}_\tau$ , which nearly coincides with the equilibrium layer of constant total stress (sum of viscous and turbulent shear stresses) as shown in figure 10(b).

In § A.2 of the Appendix, additional comparisons in the Reynolds stress components against data from Pirozzoli & Bernardini (Reference Pirozzoli and Bernardini2011), Sillero et al. (Reference Sillero, Jiménez and Moser2013) and Cogo et al. (Reference Cogo, Salvadore, Picano and Bernardini2022) at similar Reynolds numbers are included to assess the collapse of the data across different Mach number and wall-temperature conditions, as well as to demonstrate the validity of the current DNS against the established datasets in the literature (figure 37).

Figure 11. Profiles of turbulent heat-flux components $\overline {\rho T^{\prime\prime} u^{\prime\prime}_i}/\bar {\rho } T^* u^*_\tau$ as a function of Reynolds number. The shaded region indicates the region of $5\leqslant z^*\leqslant 30$ .

3.4. Turbulent heat flux and thermodynamic variables

The turbulent heat flux $\overline {\rho T^{\prime\prime}{u_i}^{\prime\prime}}$ plays an important role in Reynolds-averaged Navier–Stokes (RANS) modelling (Wilcox et al. Reference Wilcox2006; Bowersox Reference Bowersox2009; Bowersox & North Reference Bowersox and North2010). Figure 11 outlines the streamwise and wall-normal components of the turbulent heat flux in semi-local units. Here, the turbulent heat flux is normalised by the semi-local velocity $u_\tau ^*$ and the semi-local temperature $T^*$ defined as

(3.7) \begin{equation} u_\tau ^* = \sqrt {\frac {\rho _w}{\overline {\rho }}} u_\tau , \quad \quad T^* = \gamma \widetilde {T} {M^*_\tau }^2, \end{equation}

where $M_\tau ^* = u_\tau ^* / \sqrt {\gamma R {\widetilde {T}^{^{^{^{^{\!\!\!}}}}}}}$ is the semi-local friction Mach number and $\gamma$ is the specific heat ratio. The semi-local scaling leads an overall good collapse of both turbulent heat-flux components in the inner region of the boundary layer across the Reynolds number range, with $\overline {\rho T^{\prime\prime}u''}$ and $\overline {\rho T^{\prime\prime}w''}$ exhibiting a positive and negative peak at $z^* \simeq 4$ , respectively. In particular, the location of the negative peak of $\overline {\rho T^{\prime\prime}w''}$ at $z^* \simeq 4$ coincides with the positive near-wall peak location of $\partial \widetilde {T}/\partial z$ (figure 8 b), which aligns with the gradient-diffusion hypothesis commonly used for RANS modelling. When the Reynolds number ${Re}_\tau$ is increased, a noticeable increase in the magnitude of both turbulent heat-flux components occurs in the outer region of the boundary layer. However, the relative increase becomes significantly smaller when the Reynolds number ${Re}_\tau$ increases from $1947$ to $2480$ , with peak values of $\overline {\rho T^{\prime\prime}{u_i}^{\prime\prime}}/(\bar {\rho } T^* u^*_\tau )$ settling to $- 12$ at $z^* \approx 500$ and $4.6$ at $z^* \approx 2500$ for the streamwise and wall-normal components, respectively. It is worth mentioning that the streamwise turbulent heat flux also shows a local peak in the buffer layer at $z^* \approx 25$ , with the peak value converging to $\overline {\rho T^{\prime\prime}u''}/(\bar {\rho } u^*_\tau T^*) \approx -14$ at ${Re}_\tau = 2480$ , while a similar buffer layer peak is absent from the wall-normal component. The different behaviour of the streamwise and wall-normal turbulent heat fluxes in the buffer layer is likely caused by the more complex balance of the turbulent heat-flux budget terms in the near-wall region for the streamwise component than for the wall-normal component.

Figure 12. (a,b) Turbulent Prandtl number $Pr_t$ and (c,d) total Prandtl number ${\textit{Pr}}_{\textit{total}}$ as a function of Reynolds number.

The turbulent Prandtl number $Pr_t$ , which quantifies the relationship between turbulent transport of momentum and heat, is defined as

(3.8) \begin{equation} Pr_t =\frac {[\overline {\rho u''w''}](\partial \tilde {T}/\partial {z})}{[\overline {\rho T^{\prime\prime}w'' }](\partial \tilde {u}/\partial {z})}. \end{equation}

Figure 12(a) shows the turbulent Prandtl number ${\textit{Pr}}_t$ as a function of Reynolds number. Consistent with previous studies (Beresh et al. Reference Beresh, Henfling, Spillers and Pruett2011; Zhang et al. Reference Zhang, Duan and Choudhari2018; Huang et al. Reference Huang, Nicholson, Duan, Choudhari and Bowersox2020b ), $Pr_t$ has values close to unity in most of the outer region of the boundary layer. A constant value of $Pr_t = 0.85$ provides a good fit to the DNS data over $z/\delta \simeq 0.1 - 0.4$ , and this constant value is commonly assumed in gradient-diffusion models to approximate the turbulent heat-flux term within the RANS framework. For $z/\delta \simeq 0.4 - 1.0$ , the linear relation $Pr_t = 0.95-(1/4)(z/\delta )$ provides a better fit to the DNS data, with minimal Reynolds number dependence observed. Figure 12(b) shows that the assumption of a constant turbulent Prandtl number fails in the near-wall region, where $Pr_t$ undergoes an overshoot at $z^* \simeq 4$ (where the turbulent heat flux $\overline {\rho T^{\prime\prime}w'' }$ becomes zero) and shows an observable Reynolds number dependence. This trend is consistent with those reported in multiple previous studies (Huang et al. Reference Huang, Nicholson, Duan, Choudhari and Bowersox2020b , Reference Huang, Duan and Choudhari2022; Cogo et al. Reference Cogo, Baù, Chinappi, Bernardini and Picano2023). As seen in figure 11, both $\overline {\rho u'' w''}$ and $\overline {\rho T^{\prime\prime}w''}$ approach zero at the wall. This behaviour leads to a 0/0 numerical instability, which makes $Pr_t$ highly sensitive to numerical errors (Chen et al. Reference Chen, Zhu, Shi and Yang2023). Huang et al. (Reference Huang, Coleman, Spalart and Yang2023) alleviates this instability by introducing the total Prandtl number ${\textit{Pr}}_{\textit{total}}$ , defined as

(3.9) \begin{equation} {\textit{Pr}}_{\textit{total}} = \frac {1 + {\rm \mu} _t / \overline {{\rm \mu} }}{1 / {\textit{Pr}} + ({\rm \mu} _t / \overline {{\rm \mu} }) / {\textit{Pr}}_t}, \end{equation}

where the turbulent viscosity ${\rm \mu} _t$ and the turbulent Prandtl number $Pr_t$ are obtained using (3.5)and (3.8), respectively. In addition to introducing the definition of ${\textit{Pr}}_{\textit{total}}$ , Huang et al. (Reference Huang, Coleman, Spalart and Yang2023) proposed a modelled form of the total Prandtl number, $({\textit{Pr}}_{\textit{total}})_m$ , based on Reynolds-averaged momentum and energy equations and assuming constant total shear stress and heat flux in the inner regions of the TBL. Here $({\textit{Pr}}_{\textit{total}})_m$ can be written as

(3.10) \begin{equation} ({\textit{Pr}}_{\textit{total}})_m = \frac {1 + ({\rm \mu} _t / \overline {{\rm \mu} })_m}{1 / {\textit{Pr}} + ({\rm \mu} _t/\overline {{\rm \mu} })_m / ({\textit{Pr}}_t)_m}, \end{equation}

wherein the turbulent to molecular viscosity ratio ${\rm \mu} _t/\overline {{\rm \mu} }$ and $Pr_t$ are modelled as

(3.11) \begin{equation} ({\rm \mu} _t/\overline {\rm \mu} )_m = \kappa z^*\left[ 1 - \exp ( -{z^*}/{A^*})\right]^2, \qquad ({\textit{Pr}}_t)_m = 1.05-0.2\operatorname {tanh}^3({z^*}/{A^*_{\textit{Pr}}}). \end{equation}

The model constants in (3.11) were set to $A^*=17$ and $A^*_{Pr}=70$ by Huang et al. (Reference Huang, Coleman, Spalart and Yang2023) to best match the DNS results for compressible turbulent channel flows. The modelled turbulent viscosity $({\rm \mu} _t/\overline {{\rm \mu} })_m$ is compared against the DNS data in figure 10(a), while the comparison for the turbulent Prandtl number $(Pr_t)_m$ is shown in figure 12(b). Figures 12(c) and 12(d) plot the total Prandtl number versus the wall-normal distance for all Reynolds numbers. Similar to $Pr_t$ , ${\textit{Pr}}_{\textit{total}}$ is nearly a constant of $Pr_t \simeq 0.85$ in the range of $z/\delta \simeq 0.1 - 0.4$ and follows a linear relationship of ${\textit{Pr}}_{\textit{total}} = 0.95 - (1/4)(z/\delta )$ over $z/\delta \simeq 0.4 - 1.0$ for all Reynolds numbers. As expected, the definition of ${\textit{Pr}}_{\textit{total}}$ successfully accounts for the near-wall singularity and Reynolds number dependence in $Pr_t$ , leading to a constant asymptotic value of ${\textit{Pr}}_{\textit{total}} = 0.701$ for $z^* \lesssim 5$ at all Reynolds numbers. The modelling approach of (3.10) successfully predicts the general behaviour of ${\textit{Pr}}_{\textit{total}}$ in the near-wall regions for the current boundary-layer flow, which is consistent with its good performance for compressible channel flows as shown in Huang et al. (Reference Huang, Coleman, Spalart and Yang2023). Such a trend suggests the robustness of the model to flow configurations and parameters.

The fluctuating thermodynamic quantities are also relevant to turbulence modelling for high-speed flows, as they appear in advanced RANS models such as the algebraic energy flux (AEF) model by Bowersox (Reference Bowersox2009). Figure 13 plots the root-mean-square (r.m.s.) temperature, density and pressure in semi-local scaling. The temperature fluctuations $T^{\prime}_{rms}$ exhibit a double-peaked structure that is characteristic of a moderately cooled TBL (Gibis et al. Reference Gibis, Sciacovelli, Kloker and Wenzel2024). Consistent with the wall-normal turbulent heat flux $\overline {\rho T^{\prime\prime}w''}$ , we see a similarly good collapse of $\rho ^{\prime}_{rms}/\rho ^*$ and $T^{\prime}_{rms}/T^*$ in the viscous sublayer, buffer layer and at least part of the log layer at various Reynolds numbers, while the peak near the edge of the boundary layer is spread across a broad range of $z^*$ , with the outer region peak remaining to be strongly dependent on the Reynolds number. For the pressure fluctuations $p^{\prime}_{rms}$ , however, the fluctuation intensity increases with Reynolds number across the boundary layer when normalised by $\tau _w$ . Such a trend is consistent with previous DNS findings at lower Reynolds numbers (Bernardini & Pirozzoli Reference Bernardini and Pirozzoli2011; Huang et al. Reference Huang, Duan, Casper, Wagnild and Bitter2020a ; Cogo et al. Reference Cogo, Baù, Chinappi, Bernardini and Picano2023).

Figure 13. Profiles of root-mean-square thermodynamic quantities as a function of Reynolds number. Here, $\rho ^* = \gamma \overline {\rho } {M^*_\tau }^2$ is the semi-local density, and the semi-local temperature $T^*$ and friction Mach number $M_\tau ^*$ are defined in (3.7).

Figure 14. (a) Skewness, (b) flatness and (c) intermittency factor of the streamwise velocity fluctuations as a function of the Reynolds number. The diamond symbols denote the incompressible TBL DNS data of Örlü & Schlatter (Reference Örlü and Schlatter2013) at ${Re}_\tau =1280$ , and the circles correspond to the incompressible TBL data by Eitel-Amor et al. (Reference Eitel-Amor, Örlü and Schlatter2014) at ${Re}_\tau =2479$ . The horizontal line denotes Gaussian skewness and flatness values (i.e. $S(u)$ = 0 in (a) and $F(u) = 3$ in (b), respectively).

3.5. Higher-order statistics

The skewness and flatness (kurtosis) of a given flow quantity $\phi$ is defined by $S(\phi )=\overline {\phi '^3}/(\overline {\phi ^{{\prime}2}})^{3/2}$ and $F(\phi )=\overline {\phi '^4}/(\overline {\phi ^{{\prime}2}})^{2}$ , respectively, with Gaussian skewness and flatness values defined as $S(\phi ) = 0$ and $F(\phi ) = 3$ . Figures 14(a) and 14(b) plot the skewness and flatness profiles of the streamwise velocity fluctuation $u'$ as a function of Reynolds number. The DNS data of an incompressible boundary layer at ${Re}_\tau =1280$ from Örlü & Schlatter (Reference Örlü and Schlatter2013) is also included for comparison with the compressible case. The present DNS case exhibits positive skewness and flatness for $z^* \approx 1{-}15$ , after which negative skewness is observed throughout the remainder of the boundary layer. For $50 \lesssim z^* \lesssim 1000$ , $S(u) \approx 0$ and $F(u) \approx 3$ is maintained, indicating nearly Gaussian behaviour throughout a significant portion of the boundary layer. Such near-wall behaviour of the skewness and flatness of $u'$ is similar to those observed in the incompressible channel flow data by Moser, Kim & Mansour (Reference Moser, Kim and Mansour1999), which suggests the effectiveness of semi-local scaling in accounting for the Mach number and wall-cooling effects on near-wall turbulence dynamics. As the Reynolds number is increased, the sharp deviation from Gaussian values apparently becomes more localised to the edge of the boundary layer. In other words, the extent of the Gaussian region tends to become larger with increasing Reynolds number, which is consistent with the growth of the boundary layer’s logarithmic region at higher Reynolds numbers.

Previous studies of hypersonic TBLs (Duan et al. Reference Duan, Beekman and Martin2011; Huang et al. Reference Huang, Duan and Choudhari2022) found that the rapid changes in $S(u)$ and $F(u)$ in the outer edge of the boundary layer (when approaching the free stream from within the boundary layer) are indicative of the onset of intermittency that measures the wall-ward extent of incursions of the external, irrotational flow into the boundary layer. As such, the intermittency factor can be plotted in figure 14(c), defined here as the departure from the Gaussian flatness value of $3$ :

(3.12) \begin{equation} \gamma = \frac {3}{F(u)}. \end{equation}

The intermittency factor tends to remain at or around $\gamma =1$ for a longer wall-normal extent for the higher-Reynolds-number profiles, indicating a delay in the onset of intermittency with Reynolds number. Interestingly, the Mach number seems to have caused an even larger delay in the onset of intermittency, as shown when comparing the incompressible data of Eitel-Amor, Örlü & Schlatter (Reference Eitel-Amor, Örlü and Schlatter2014) with the current hypersonic data at a common Reynolds number of ${Re}_\tau \approx 2500$ .

3.6. Budget terms

In this section the budget terms in the transport equations for TKE $\overline {\rho }\widetilde {k_t} = \overline {\rho u^{\prime\prime}_iu^{\prime\prime}_i}/2$ (3.13), temperature variance $\overline {\rho T^{{\prime\prime}2}}$ (3.15) and turbulent heat flux $\overline {\rho T^{\prime\prime}{u_i}^{\prime\prime}}$ (3.18) are decomposed into their respective budget terms to better understand the turbulent dynamics.

3.6.1. Turbulent kinetic energy and temperature variance budget

The budget for TKE $\overline {\rho }\widetilde {k_t} = \overline {\rho u^{\prime\prime}_iu^{\prime\prime}_i}/2$ is comprised of various source and sink terms that balance the following transport equation according to the formulation by Wilcox et al. (Reference Wilcox2006):

(3.13) \begin{equation} \overline {\rho }\frac {D \widetilde {k_t}}{D t}=\frac {\partial \left(\overline {\rho }\widetilde {k_t}\right)}{\partial {t}}+\frac {\partial }{\partial {x_k}} \left(\overline {\rho }\widetilde {u}_k\widetilde {k_t} \right)=P_{k_t}+T_{k_t}+\Pi _{k_t}-\overline {\rho }\epsilon _{k_t}+D_{k_t}+M_{k_t} \end{equation}

with

(3.14a) \begin{align} P_{k_t} &= -\overline {\rho u^{\prime\prime}_iu^{\prime\prime}_k}\frac {\partial \widetilde {u}_{i}}{\partial x_k} ,\end{align}

(3.14b) \begin{align} T_{k_t} &=-\frac {\partial }{\partial x_k}\left (\frac {1}{2}\overline {\rho u^{\prime\prime}_iu^{\prime\prime}_iu^{\prime\prime}_k} \right ), \end{align}

(3.14c) \begin{align} \Pi _{k_t} &= \pi _{k_t}^t+\pi _{k_t}^d = -\frac {\partial }{\partial x_i}\left(\overline {p'u^{\prime\prime}_i}\right)+\overline {p'\frac {\partial u^{\prime\prime}_i}{\partial x_i}} ,\end{align}

(3.14d) \begin{align} \overline {\rho }\epsilon _{k_t} &= \overline {t^{\prime}_{ik}\frac {\partial u^{\prime\prime}_i}{\partial x_k}} ,\end{align}

(3.14e) \begin{align} D_{k_t} &= \frac {\partial }{\partial x_k}\left(\overline {t^{\prime}_{ik}u^{\prime\prime}_i}\right), \end{align}

(3.14f) \begin{align} M_{k_t} &=\overline {u^{\prime\prime}_i}\left (\frac {\partial \overline {t_{ik}}}{\partial x_k}-\frac {\partial \overline {p}}{\partial x_i} \right ), \end{align}

where $P_{k_t}$ is the production term, $T_{k_t}$ is the turbulent transport term, $\Pi _{k_t}$ is the pressure term (pressure diffusion and pressure dilatation), $-\overline {\rho }\epsilon _{k_t}$ is the viscous dissipation per unit volume, $D_{k_t}$ is the viscous diffusion and $M_{k_t}$ is the mass flux contribution associated with density fluctuations. Here, $t_{ik}=2{\rm \mu} [S_{ik}-(1/3)S_{qq}\delta _{ik}]$ is the viscous stress tensor and $S_{ik}=0.5(\partial u_i/\partial x_k + \partial u_k/\partial x_i)$ is the rate of strain tensor.

Figure 15 plots each term of the TKE budget in semi-local units. The compressible DNS of Cogo et al. (Reference Cogo, Salvadore, Picano and Bernardini2022) is included for comparison purposes. In the near-wall region ( $z^* \lesssim 5$ ), TKE transport is mainly governed by a balance between viscous diffusion $D_{k_t}$ and viscous dissipation $-\overline {\rho } \epsilon _{k_t}$ . Further from the wall ( $5 \lesssim z^* \lesssim 70$ ), TKE transport is mainly dominated by the production term $P_{k_t}$ , with loss contributions coming from viscous diffusion $D_{k_t}$ , turbulent transport $T_{k_t}$ and viscous dissipation $-\overline {\rho } \epsilon _{k_t}$ . The mass flux term $M_{k_t}$ and the pressure diffusion/dilatation term $\Pi _{k_t}$ minimally contribute to the overall transport of TKE throughout the entire boundary layer. As shown, the Reynolds number has minimal effects on the TKE dynamics, and discrepancies from the data of Cogo et al. (Reference Cogo, Salvadore, Picano and Bernardini2022) are nearly indistinguishable when the semi-local scaling is used.

Figure 15. Budgets of the transport equation for TKE: (a) production term $P_{k_t}$ , (b) turbulent transport $T_{k_t}$ , (c) pressure diffusion/pressure dilatation term $\Pi _{k_t}$ , (d) viscous dissipation per unit volume $-\overline {\rho } \epsilon _{k_t}$ , (e) viscous diffusion $D_{k_t}$ , ( f) mass flux contribution due to density fluctuations $M_{k_t}$ . Variables are normalised by $\overline {\rho }u^{*3}_{\tau }/z^*_\tau$ and plotted in the semi-local scale $z^*$ .

The transport equation for the temperature variance $\overline {\rho T^{{\prime\prime}2}}$ is further analysed to shed light on the underlying mechanisms that govern the thermal transport. The transport equation for $\overline {\rho T^{{\prime\prime}2}}$ is written as follows (Gerolymos & Vallet Reference Gerolymos and Vallet2014):

(3.15) \begin{equation} \frac {D\overline {\rho T^{{\prime\prime}2}}}{Dt}=\frac {\partial \overline {\rho T^{{\prime\prime}2}}}{\partial t}+\frac {\partial \overline {\rho T^{{\prime\prime}2}}\widetilde {u}_k}{\partial x_k}=D_{T}+P_{T}-\overline {\rho }\epsilon _{T}+C_{T}+T_{T}+B_{T} \end{equation}

with

(3.16a) \begin{align} D_{T} &= -\frac {\partial \overline {\rho {T^{\prime\prime}}^2 {u_k}^{\prime\prime}}}{\partial x_k}-\frac {2}{c_v}\frac {\partial \overline {T'{q_k}^{\prime}}}{\partial x_k}, \end{align}

(3.16b) \begin{align} P_{T} &= -2\overline {\rho T^{\prime\prime}{u_k}^{\prime\prime}}\frac {\partial \widetilde {T}}{\partial x_k}+\frac {2}{c_v}\overline {T't_{jk}^{\prime}}\overline {S}_{jk}+\frac {2}{c_v}\overline {T'S_{jk}^{\prime}}\overline {t}_{jk}-\frac {2}{c_v}\overline {T'p'}\frac {\partial \widetilde {u}_k}{\partial x_k} ,\end{align}

(3.16c) \begin{align} \overline {\rho }\epsilon _{T} &= -\frac {2}{c_v}\overline {q_k'\frac {\partial T'}{\partial x_k}}, \end{align}

(3.16d) \begin{align} C_{T} &= \frac {2}{c_v}\overline {T^{\prime\prime}}\left (\overline {t_{jk}S_{jk}}-\overline {p\frac {\partial u_k}{\partial x_k}}-\frac {\partial q_k}{\partial x_k}\right ), \end{align}

(3.16e) \begin{align} T_{T} &= \frac {2}{c_v}\overline {T't^{\prime}_{jk}S^{\prime}_{jk}} ,\end{align}

(3.16f) \begin{align} B_{T} &= -\frac {2}{c_v}\overline {pT'\frac {\partial u^{\prime\prime}_k}{\partial x_k}}, \end{align}

where $S_{ij}$ is the strain-rate tensor, $q_i=-\kappa ({\partial T}/{\partial x_i})$ is the molecular heat flux, $\kappa$ is the heat conductivity, $c_v$ is the heat capacity at constant volume, $D_{T}$ is the diffusion, $P_{T}$ is the production, $-\overline {\rho }\epsilon _{T}$ is the dissipation, $C_{T}$ is the contribution of compressibility, $T_{T}$ is the contribution of triple correlation and $B_{T}$ is the contribution of fluctuating dilatation $ {\partial u^{\prime\prime}_k}/{\partial x_k}$ .

Figure 16. Budgets of the transport equation for the temperature variance $\overline {\rho T^{{\prime\prime}2}}$ : (a) the diffusion term $D_{T}$ , (b) the production term $P_{T}$ , (c) the dissipation term $-\overline {\rho }\epsilon _{T}$ , (d) the compressibility terms $C_{T}$ , (e) the triple correlation term $T_{T}$ , ( f) the fluctuating dilatation term $B_{T}$ . Variables are normalised by $\overline {\rho } u_\tau ^* \widetilde {T}^2/z_\tau ^*$ and plotted in the semi-local scale $z^*$ .

Figure 17. Budget terms of temperature variance $\overline {\rho T^{{\prime\prime}2}}$ at ${Re}_\tau = 2480$ . Variables are normalised by $\overline {\rho } u_\tau ^* \widetilde {T}^2/z_\tau ^*$ and plotted in the semi-local scale $z^*$ . The r.m.s. of temperature fluctuations is plotted on a separate $y$ axis to indicate the local peak location. The shaded region indicates the region of $5\leqslant z^*\leqslant 30$ .

Figure 16 presents the budgets for the temperature variance $\overline {\rho T^{{\prime\prime}2}}$ as a function of Reynolds number. Unlike the TKE budget terms that demonstrate a consistent collapse across different Reynolds numbers when plotted against $z^*$ (figure 15), the temperature variance budgets show a noticeable decrease in magnitude when the Reynolds number ${Re}_\tau$ is increased. However, the relative decrease becomes significantly smaller when the Reynolds number ${Re}_\tau$ increases from $1947$ to $2480$ . Figure 17 further shows the temperature variance budgets at ${Re}_\tau = 2480$ along with the r.m.s. temperature $T^{\prime}_{rms}/T^*$ . In the near-wall region ( $z^* \lesssim 40$ ), the dynamics of the temperature variance is dominated by a balance among diffusion $D_{T}$ , production $P_{T}$ , dissipation $-\overline {\rho }\epsilon _{T}$ and fluctuating dilatation $B_T$ . This region also features a sharp rise in the r.m.s. temperature fluctuations, resulting in a pronounced local peak at $z^* \simeq 30$ across all Reynolds numbers. In this context, the production term $P_T$ , which represents the conversion of mean thermal energy into thermal fluctuations, plays a central role. The alignment of the localised peak location in $P_T$ and r.m.s. temperature fluctuations indicates that increased production amplifies temperature variance. A similar observation was reported by Cogo et al. (Reference Cogo, Baù, Chinappi, Bernardini and Picano2023), who also found that the near-wall peak in temperature variance coincides with enhanced thermal production, reinforcing the connection between production mechanisms and intensified temperature fluctuations in this region.

While the budget analysis highlights the key mechanisms governing the transport of TKE and temperature variance, RANS models such as the AEF model of Bowersox (Reference Bowersox2009) depend on the characteristic turbulent velocity and thermal time scales to approximate eddy viscosity and diffusivity without solving the full transport equations. These time scales are particularly crucial for capturing the turbulent dynamics of compressible boundary layers subjected to wall cooling. Specifically, Bowersox (Reference Bowersox2009) defined the following characteristic turbulent velocity scale $\tau _u$ and thermal time scale $\tau _T^B$ in his derivation of the AEF model:

(3.17) \begin{eqnarray} \tau _u = \left(\overline {\rho }\widetilde {k_t}\right) / \left(\overline {\rho }\epsilon _{k_t}\right)\!,\!\!\!\! \quad \quad \tau _T^B = 2 \left [(P_{k_t} - \overline {\rho }\epsilon _{k_t})/\overline {\rho }\widetilde {k}_t + (P_T + \overline {\rho }\epsilon _T )/ \frac {1}{2}\overline {\rho }\widetilde {T^{{\prime\prime}2}} \right ]^{-1}\!, \end{eqnarray}

where $P_{k_t}$ and $\overline {\rho }\epsilon _{k_t}$ are production and dissipation terms, respectively, in the budget for TKE $\overline {\rho }\widetilde {k_t} = \overline {\rho u^{\prime\prime}_iu^{\prime\prime}_i}/2$ (3.13), and $P_T$ and $\overline {\rho }\epsilon _T$ are production and dissipation terms, respectively, in the budget for the temperature variance $\overline {\rho T^{{\prime\prime}2}}$ (3.15). Bowersox (Reference Bowersox2009) further introduced a model constant $\sigma _\nu \equiv \tau ^B_T/\tau _u \approx 0.2$ by assuming the thermal time scale $\tau ^B_T$ to be of the same order as the velocity time scale $\tau _u$ across the boundary layer.

To validate or improve the modelling assumptions of the AEF model, figure 18(a) shows the turbulent thermal time scale $\tau _T^B$ , which represents the persistence of thermal fluctuations before they dissipate. This time scale increases with $z^*$ , reflecting reduced dissipation away from the wall. Similarly, figure 18(b) presents the velocity time scale, $\tau _u$ , which also increases with $z^*$ , though the growth rate differs, particularly in the buffer and log layers. Figure 18(c) plots the ratio of the thermal time scale to the velocity time scale $\sigma _\nu$ , which shows that $\sigma _\nu$ varies significantly with the wall-normal distance, specifically in the near-wall region where dissipation dominates turbulent transport. However, in the logarithmic region ( $50\lesssim z^* \lesssim 350$ ) and in the outer layer, the constant value of $\sigma _\nu$ approaches approximately $0.2$ . This behaviour suggests that assuming a constant $\sigma _\nu$ throughout the boundary layer may be inappropriate, particularly in near-wall regions. Taking into account the spatial variation in $\tau _T^B$ and $\tau _u$ may therefore be essential to accurately model turbulent heat flux and eddy diffusivity in compressible boundary layers.

Figure 18. (a) Turbulent thermal time scale $\tau _T^B = 2 [(P_{k_t} - \overline {\rho }\epsilon _{k_t})/\overline {\rho }\widetilde {k}_t + (P_T + \overline {\rho }\epsilon _T )/ ({1}/{2})\overline {\rho }\widetilde {T^{{\prime\prime}2}} ]^{-1}$ , (b) velocity time scale $\tau _u = \overline {\rho }\widetilde {k_t} / \overline {\rho }\epsilon _{k_t}$ and (c) model constant $\sigma _\nu =\tau _T^B/\tau _u$ as a function of wall-normal distance $z^*$ .

3.6.2. Turbulent heat-flux budget

To gain a detailed insight into the underlying mechanisms governing the turbulent heat flux, we analyse here the decomposition of $\overline {\rho T^{\prime\prime}{u_i}^{\prime\prime}}$ into its respective budget terms. Furthermore, the budget terms are particularly relevant to the formulation of the AEF model developed by Bowersox (Reference Bowersox2009) and Bowersox & North (Reference Bowersox and North2010). Existing literature suggests that the extension of this model to hypersonic, cold-wall TBLs is promising, and the analysis of the heat-flux budget terms using high-fidelity datasets is paramount to confirming the model assumptions (Nicholson et al. Reference Nicholson, Huang, Duan, Choudhari, Morreale and Bowersox2022). The formulation of the budget terms is outlined here as (Shahab et al. Reference Shahab, Lehnasch, Gatski and Comte2011; Nicholson et al. Reference Nicholson, Huang, Duan, Choudhari, Morreale and Bowersox2022)

(3.18) \begin{equation} \frac {D\overline {\rho T^{\prime\prime}{u_i}^{\prime\prime}}}{Dt}=\frac {\partial \overline {\rho T^{\prime\prime}{u_i}^{\prime\prime}}}{\partial t}+\frac {\partial \overline {\rho T^{\prime\prime}{u_i}^{\prime\prime}}\widetilde {u}_k}{\partial x_k}=P_{Tu_i}^1+P_{Tu_i}^2+\Pi _{Tu_i}-\overline {\rho }\epsilon _{Tu_i}+D_{Tu_i}+C_{Tu_i}, \end{equation}

with

(3.19a) \begin{align} P_{Tu_i}^1 &= -\overline {\rho T^{\prime\prime}{u_k}^{\prime\prime}}\frac {\partial \widetilde {u}_i}{\partial x_k}-\overline {\rho {u_i}^{\prime\prime}{u_k}^{\prime\prime}}\frac {\partial \widetilde {T}}{\partial x_k}+\frac {1}{c_v}\overline {{u_i}^{\prime\prime}t_{jk}^{\prime}}\widetilde {S}_{kj}-(\gamma -1)\overline {\rho T^{\prime\prime}{u_i}^{\prime\prime}}\widetilde {S}_{kk}, \end{align}

(3.19b) \begin{align} P_{Tu_i}^2 &= \frac {1}{c_v}\left (\overline {{u_i}^{\prime\prime}t_{jk}^{\prime}S_{kj}^{\prime\prime}}-\overline {{u_i}^{\prime\prime}p'S_{kk}^{\prime\prime}}\right ), \end{align}

(3.19c) \begin{align} \Pi _{Tu_i}&=-\overline {T^{\prime\prime}\frac {\partial p}{\partial x_i}}, \end{align}

(3.19d) \begin{align} \overline {\rho }\epsilon _{Tu_i} &= -\frac {1}{c_v}\overline {{q_k}^{\prime}\frac {\partial {u_i}^{\prime\prime}}{\partial x_k}}+\overline {t_{ik}^{\prime}\frac {\partial T^{\prime\prime}}{\partial x_k}} ,\end{align}

(3.19e) \begin{align} D_{Tu_i} &= -\frac {\partial \overline {\rho T^{\prime\prime}{u_i}^{\prime\prime}{u_k}^{\prime\prime}}}{\partial x_k}+\frac {\partial \overline {T^{\prime\prime}t_{ik}^{\prime}}}{\partial x_k}-\frac {1}{c_v}\frac {\partial \overline {{u_i}^{\prime\prime}{q_k}^{\prime}}}{\partial x_k}, \end{align}

(3.19f) \begin{align} C_{Tu_i} &= \overline {T^{\prime\prime}}\frac {\partial \overline {t}_{ik}}{\partial x_k}+\frac {1}{c_v} \left[\overline {{u_i}^{\prime\prime}}\overline {t}_{jk}\widetilde {S}_{kj}+\overline {{u_i}^{\prime\prime}S_{kj}^{\prime\prime}}(\overline {t}_{kj}-\overline {p}\delta _{kj}) \right]. \end{align}

Here, $P_{Tu_i}^1$ is the production due to the mean temperature and velocity gradients, $P_{Tu_i}^2$ is the production due to the fluctuating strain rate, $\Pi _{Tu_i}$ is the ‘pressure-scrambling’ term, $-\overline {\rho }\epsilon _{Tu_i}$ is the turbulent viscous-thermal dissipation term, $D_{Tu_i}$ is the turbulent viscous-thermal transport contribution and $C_{Tu_i}$ is the term associated with compressibility.

Figure 19. Streamwise turbulent heat-flux budget terms: (a) the first production term $P_{\textit{Tu}}^1$ , (b) the second production term $P_{\textit{Tu}}^2$ , (c) the turbulent viscous-thermal dissipation term $-\overline {\rho }\epsilon _{\textit{Tu}}$ , (d) the compressibility terms $C_{\textit{Tu}}$ , (e) the ‘pressure-scrambling term’ $\Pi _{\textit{Tu}}$ and ( f) the turbulent-viscous transport term $D_{\textit{Tu}}$ . Variables are normalised by $\overline {\rho }u_\tau ^{*2}\widetilde {T}/z_\tau ^*$ and plotted in semi-local scale $z^*$ .

Figure 19 outlines the budget terms for the streamwise component of the turbulent heat flux $\overline {\rho T^{\prime\prime}u''}$ as a function of Reynolds number. Very near the wall ( $z^* \lesssim 2$ ), the viscous-thermal transport term $D_{\textit{Tu}}$ acts as a sink for the turbulent heat flux, while the viscous-thermal dissipation term $\rho \epsilon _{\textit{Tu}}$ counteracts in the form of an equal but opposite source. This finding is contrary to a Mach 2 adiabatic wall case reported by Shahab et al. (Reference Shahab, Lehnasch, Gatski and Comte2011), which displays zero contributions from these terms, indicating a clear sensitivity to wall cooling. Away from the wall in the buffer layer ( $5 \lesssim z^* \lesssim 15$ ), the first production term ( $P_{\textit{Tu}}^1$ ) increases from zero to slightly positive, then drops to a stronger negative value. Terms such as the viscous-thermal dissipation term $\overline {\rho } \epsilon _{\textit{Tu}}$ , the viscous-thermal transport term $D_{\textit{Tu}}$ and the compressibility term $C_{\textit{Tu}}$ act in a neutral balance with the production, depending on the location throughout the boundary layer. Lastly, the second production terms ( $P_{\textit{Tu}}^2$ ) and the pressure-scrambling term ( $\Pi _{\textit{Tu}}$ ) have very small contributions throughout the entire boundary layer.

Figure 20. Wall-normal turbulent heat-flux budget terms: (a) the first production term $P_{Tw}^1$ , (b) the second production term $P_{Tw}^2$ , (c) the turbulent viscous-thermal dissipation term $-\overline {\rho }\epsilon _{\textit{Tw}}$ , (d) the compressibility terms $C_{Tw}$ , (e) the ‘pressure-scrambling term’ $\Pi _{Tw}$ , and ( f) the turbulent-viscous transport term $D_{Tw}$ . Variables are normalised by $\overline {\rho }u_\tau ^{*2}\widetilde {T}/z_\tau ^*$ and plotted in semi-local scale $z^*$ .

Figure 20 further outlines the budget terms for the wall-normal component of the turbulent heat flux as a function of Reynolds number. The dynamics of the wall-normal heat-flux transport $\overline {\rho T^{\prime\prime}w''}$ are largely governed by a balance between the first production term ( $P_{Tw}^1$ ) and the pressure-scrambling term ( $\Pi _{Tw}$ ) in the buffer and log layers, while the relative contributions from the other terms are small. Interestingly, the viscous-thermal transport term ( $D_{Tw}$ ) and the viscous-thermal dissipation term ( $\overline {\rho } \epsilon _{Tw}$ ) do not manifest in the wall-normal heat-flux component very near the wall ( $z^* \lesssim 2$ ), which follows a drastically different behaviour from those in the streamwise components as shown in figures 19(c) and 19( f). The relative contribution from the compressibility term $C_{Tw}$ is practically zero throughout the entire boundary layer, which is unlike the compressibility term in the streamwise heat-flux component ( $C_{\textit{Tu}}$ ) that shows a non-negligible contribution in the buffer region (figure 19 d). The different balances of the streamwise and wall-normal heat-flux budget terms will likely lead to the different near-wall behaviours of the streamwise and wall-normal heat-flux components as shown in figure 11.

Similar to the budget terms for the temperature variance $\overline {\rho T^{{\prime\prime}2}}$ , the budgets for $\overline {\rho T^{\prime\prime}u''}$ and $\overline {\rho T^{\prime\prime}w''}$ generally show a noticeable reduction in peak levels when the Reynolds number increases from ${Re}_\tau = 1116$ to ${Re}_\tau = 1947$ . However, the clear difference becomes significantly smaller when ${Re}_\tau$ is further increased to $2480$ . This Reynolds number dependence may have led to a similar change in the turbulent heat flux as shown in figure 11.

Figure 21. Turbulent heat-flux budget terms at ${Re}_\tau = 2480$ , along with the mean temperature $\tilde {T}/T_w$ and turbulent heat-flux components $\overline {\rho T^{\prime\prime} {u_i}^{\prime\prime}}/\bar {\rho } T^* u^*_\tau$ . The budget terms are normalised by $\overline {\rho }u_\tau ^{*2}\widetilde {T}/z_\tau ^*$ and plotted in semi-local scale $z^*$ . The mean temperature and turbulent heat-flux components are plotted on separate $y$ axes. The shaded region indicates the region of $5\leqslant z^*\leqslant 30$ .

Figure 21 plots the dominant heat-flux budget terms at ${Re}_\tau = 2480$ along with the mean temperature $\widetilde {T}/T_w$ and the turbulent heat-flux components $\overline {\rho T^{\prime\prime} u^{\prime\prime}_i}/\bar {\rho } T^* u^*_\tau$ . For the streamwise component (figure 21 a), the peak of $P_{\textit{Tu}}^1$ at $z^* \simeq 4$ closely aligns with the maximum of the streamwise heat flux $\overline {\rho T^{\prime\prime}u''}$ , indicating a strong link between local production and streamwise heat flux in this region. At this same location, the diffusion term $D_{\textit{Tu}}$ exhibits a local maximum, and the compressibility term $C_{\textit{Tu}}$ also reaches its peak, highlighting the combined influence of multiple mechanisms. This alignment of production, diffusion and compressibility peaks marks a key region of energy transfer and sheds light on the enhanced near-wall heat transport under wall-cooled conditions. The dominant budget terms $D_{Tw}$ and $\overline {\rho } \epsilon _{Tw}$ for the wall-normal component exhibit peak at locations that align with the peak location of the mean temperature $\widetilde {T}/T_w$ in the buffer region (figure 21 b), highlighting this area as the key zone for the dynamics of wall-normal turbulent heat flux and turbulent thermal transport.

Figure 22 further compares the relative dominance of the four production mechanisms in $P_{Tu_i}^1$ , including those associated with the mean velocity gradient ( $-\overline {\rho T^{\prime\prime}{u_k}^{\prime\prime}}({\partial \widetilde {u}_i}/{\partial x_k})$ ), mean temperature gradient ( $-\overline {\rho {u_i}^{\prime\prime}{u_k}^{\prime\prime}}({\partial \widetilde {T}}/{\partial x_k})$ ), mean strain rate ( $({1}/{c_v})\overline {{u_i}^{\prime\prime}t_{jk}^{\prime}}\widetilde {S}_{kj}$ ) and mean velocity divergence ( $-(\gamma -1)\overline {\rho T^{\prime\prime}{u_i}^{\prime\prime}}\widetilde {S}_{kk}$ ). For the streamwise component of the turbulent heat flux, the contributions from the terms associated with the mean velocity gradient, mean temperature gradient and mean strain rate are significant. In contrast, the term linked to the mean velocity divergence is negligible. Among these mechanisms, the production term arising from the mean strain rate is the most dominant. However, this term has been omitted in the derivation of the AEF model (Bowersox Reference Bowersox2009; Bowersox & North Reference Bowersox and North2010). While this omission is reasonable for TBLs at moderate Mach numbers with adiabatic wall conditions, it may lead to inaccuracies near the wall under high-Mach-number, cold-wall conditions, as shown in figure 22. For the wall-normal heat flux, our analysis reveals that the dominant contribution to the total production originates from the mean temperature gradient, while the other three production mechanisms have negligible effects on this component. This highlights the varying significance of individual production terms depending on the heat-flux component and flow conditions.

The dynamics shown in the above analysis for both the streamwise and wall-normal turbulent heat-flux components demonstrates some of the challenges facing turbulence modelling. While the wall-normal component displays a slightly less complicated balance than the streamwise component, the underlying disparity between them makes vector-tensor based representations of the turbulent heat flux difficult (Shahab et al. Reference Shahab, Lehnasch, Gatski and Comte2011). A higher-order model that directly solves the corresponding Reynolds stress transport equations could be used. However, the model would have to be carefully tuned to demonstrate the complex interplay displayed above.

Figure 22. Decomposition of the production term of the (a) streamwise and (b) wall-normal turbulent heat-flux budget at ${Re}_\tau =2480$ . Variables are normalised by $\overline {\rho } {u_\tau ^*}^2 \widetilde {T} /z^*_\tau$ and plotted in the semi-local scale $z^*$ .