2. Observation: abatement of turbulent intensities and VLSMs

As a preliminary impression of the VLSMs, we first perform implicit large eddy simulation (LES) of turbulent boundary layers over adiabatic walls at the friction Reynolds number $Re_\tau \approx 1000$ and the free stream Mach number $M_\infty$ ranging from $2.0$ to $5.0$, governed by the Navier–Stokes equations of perfect Newtonian gases. Hereinafter, the velocity in the $x$ (streamwise), $y$ (wall-normal) and $z$ (spanwise) directions are represented by $u$, $v$ and $w$, respectively. The density, pressure and temperature are denoted by $\rho$, $p$ and $T$, related by the state equation of perfect gases. The viscous stresses and the heat transfer are obtained by the constitutive equations and Fourier's law of Newtonian fluids, with $\mu$ and $\kappa$ being the viscosity and heat conductivity.

The simulation is performed in a rectangular box with the sizes of $(80,9,10)\delta _{in}$ in $x$, $y$ and $z$ directions, discretized by $(2400, 330, 600)$ grids in the three directions, with $\delta _{in}$ the nominal boundary layer thickness at the inlet of the computational domain. At the flow inlet, the synthetic turbulence is given by the recycling-rescaling method at the reference station $x=53 \delta _{in}$ (Pirozzoli, Bernardini & Grasso Reference Pirozzoli, Bernardini and Grasso2010). The no-slip and no-penetration for velocity and isothermal condition for temperature $T_w=T_r$ ($T_r$ the recovery temperature at the given Mach number) are applied at the wall. The non-reflection conditions are enforced at the upper and outlet boundaries to eliminate the possible numerical errors that could contaminate the flow field within the computational domain. Periodic conditions are adopted in the spanwise direction.

The flow parameters are listed in table 1. The simulations are performed using the open-source code developed by Bernardini et al. (Reference Bernardini, Modesti, Salvadore and Pirozzoli2021). The convective terms are approximated by the hybrid scheme, namely the low-dissipative sixth-order kinetic preserving scheme (Pirozzoli Reference Pirozzoli2010) in the smooth region and the fifth-order weighted essentially non-oscillatory (WENO) scheme (Shu & Osher Reference Shu and Osher1988) at the flow discontinuity detected by Ducro's sensor (Ducros et al. Reference Ducros, Ferrand, Nicoud, Weber, Darracq, Gacherieu and Poinsot1999). The viscous terms are expanded into the Laplacian form and then approximated by the sixth-order central scheme. The time advancement is achieved by the third-order low-storage Runge–Kutta scheme (Wray Reference Wray1990). The subgrid stresses and heat fluxes are incorporated by the overly dissipative WENO scheme, thus no subgrid models are implemented (Kokkinakis & Drikakis Reference Kokkinakis and Drikakis2015; Ritos et al. Reference Ritos, Kokkinakis, Drikakis and Spottswood2017). The simulations have been run for more than $500 \delta /U_\infty$ before the collection of turbulent statistics, and another period of approximately $800 \delta /U_\infty$ in each case for investigation. Herein, the boundary layer thickness $\delta$ is the off-wall distance where the mean velocity reaches $99\,\%$ of the free stream values. At the streamwise station with the friction Reynolds number $Re_\tau$ of approximately $1000$, the grid resolution in the streamwise and spanwise directions are $\Delta x^+ \approx 20$ and $\Delta z^+ \approx 10$, and the first grid distance off the wall being $\Delta y^+_w \approx 0.7$ (refer to table 1), satisfying the requirement of LES. Such a method has been proven to be capable of providing acceptable results under the presently considered $Re_\tau$, grid intervals and numerical schemes (De Vanna et al. Reference De Vanna, Baldan, Picano and Benini2023).

Table 1. Flow parameters. Here, $M_\infty$ is the free stream Mach number, $Re_\infty$, $Re_\theta$ and $Re_\tau$ are the free stream, momentum and friction Mach number at the reference station. Additionally, $T_w$ is the wall temperature and $\Delta x_i$ is the mesh interval along the $x_i$ coordinate.

In figure 1(a), we present the wall-normal distribution of the van Driest transformed mean velocity $\tilde u^+_{VD}$, with the symbol $\tilde {\cdot }$ denoting its density-weighted average, as well as those by Bernardini & Pirozzoli (Reference Bernardini and Pirozzoli2011) at $M_\infty = 2$ and $Re_\tau \approx 1000$ for comparison. The results of case M2 are generally consistent with the reference data across the boundary layer, suggesting the validity of the present numerical methods. For all the cases considered, the mean velocity obeys the linear law in the viscous sublayer and the logarithmic law within $y^+=30$ and $y=0.3\delta$, but are lowered in the wake region with the rising Mach number. The density-weighted root-mean-square (r.m.s.) of velocity fluctuations $u''^*_{rms}$ shown in figure 1(b) (with ${\cdot }''$ marking the fluctuations and the superscripts $*$ the non-dimensionalization as $u''_{rms}/\sqrt {\tau _w /\rho _w}$) also manifest considerable agreement with the reference. The $u''^*_{rms}$ reaches the maximum at $y^+ \approx 12$ in the buffer layer, beyond which it gradually decays. Within a small range of $y^+ = 100 \sim 300$, the $u''^*_{rms}$ decays logarithmically in each case (Marusic et al. Reference Marusic, Monty, Hultmark and Smits2013; Lozano-Durán & Jiménez Reference Lozano-Durán and Jiménez2014), but the slopes and intercepts of the logarithmic laws (if any) are different. There exhibits a monotonic abatement as the Mach number increases.

Figure 1. Wall-normal distribution of (a) van Driest transformed mean velocity $\tilde {u}^+_{VD}$ and (b) density-weighted r.m.s. of velocity fluctuations $u''^*_{rms}$. Reference: Bernardini & Pirozzoli (Reference Bernardini and Pirozzoli2011) at $M_\infty =2$ and $Re_\tau =1000$.

It is expected that the decrease of the fluctuation intensities with the Mach numbers in the outer layer above can be associated with the flow structures. As shown in figure 2(a), the low-speed regions in the instantaneous field reach higher wall-normal locations. The spanwise length scale of these velocity streaks within the buffer layer is estimated to be $100 \delta _\nu$ ($\delta _\nu$ is the viscous length scale) (Flores & Jiménez Reference Flores and Jiménez2010; Hwang Reference Hwang2013; Yu, Xu & Pirozzoli Reference Yu, Xu and Pirozzoli2019; Cogo et al. Reference Cogo, Salvadore, Picano and Bernardini2022) and that of the low momentum region in the outer region, known as the VLSMs, is approximately the boundary layer thickness $\delta$ (Hutchins & Marusic Reference Hutchins and Marusic2007; Jiménez Reference Jiménez2018; Yu & Xu Reference Yu and Xu2022; Yu et al. Reference Yu, Fu, Tang, Yuan and Xu2023). These energetic structures are manifested as the two peaks in the pre-multiplied spanwise spectra in figure 2(b) (Huang et al. Reference Huang, Duan and Choudhari2022; Yu & Xu Reference Yu and Xu2022). Qualitatively, the distributions of the spectra and instantaneous flow structures are marginally dependent on the Mach number. However, the specific values of the inner and outer peaks extracted from each spectra, as shown in figure 2(c), suggest that the velocity streaks in the inner region are slightly strengthened but the VLSMs in the outer regions are significantly weakened with the increment of the Mach number, as also indicated by the lower spectra intensity in the outer region shown in figure 2(b), consistent with the results reported by Cogo et al. (Reference Cogo, Salvadore, Picano and Bernardini2022) and Huang et al. (Reference Huang, Duan and Choudhari2022).

Figure 2. (a) Instantaneous streamwise velocity $u/U_\infty$ in case M3, (b) the density-weighted spanwise spectra of case M3 (flooded and black solid lines) and of case M5 (red dashed lines), (c) the values of the inner and outer peaks at various Mach numbers.

Integrating the spanwise spectra gives the r.m.s. of the velocity fluctuations. The results are reported in figure 3, in which the large- and small-scale portions are separated by $\lambda _z = 0.5\delta$. Expectedly, the density-weighted r.m.s. of the small-scale velocity $u''$ is weakly dependent on the Mach number. They are also well-collapsed with the low-Reynolds-number statistics ($M_\infty = 2$, $Re_\tau = 250$) reported by Bernardini & Pirozzoli (Reference Bernardini and Pirozzoli2011) within $y^+ \approx 100$ where the VLSMs are absent, suggesting that the small-scale motions in the inner region are Mach number independent at the presently considered moderate Reynolds number and are statistically equivalent to the low-Reynolds-number turbulence. This is consistent with our previous findings with the aid of minimal flow units in turbulent channels (Yin, Huang & Xu Reference Yin, Huang and Xu2017; Yu et al. Reference Yu, Fu, Tang, Yuan and Xu2023). The large-scale portion, however, manifests consistent decrement in magnitude with the Mach number. Based on the statistics above, the abatement of the velocity fluctuation $u''^*_{rms}$ in the outer region should be ascribed to the weaker VLSMs in higher Mach number flows at the similar friction Reynolds number $Re_\tau$. Notably, when normalized by their peak values at $y=0.2\delta$, the profiles can be collapsed to a single curve in the inner region (inset of figure 3b), indicating that the superposition of the large-scale motions on the inner region (Mathis et al. Reference Mathis, Hutchins and Marusic2009; Mathis, Hutchins & Marusic Reference Mathis, Hutchins and Marusic2011) remains unaltered quantitatively.

Figure 3. Density-weighted r.m.s. of $u''$: (a) small-scale portion $u''^*_{s,rms}$; (b) large-scale portion $u''^*_{l,rms}$, where the inset shows values normalized at $y=0.2\delta$. Symbols in panel (a): Bernardini & Pirozzoli (Reference Bernardini and Pirozzoli2011) at $M_\infty =2$ and $Re_\tau =250$.

3. Dynamic relevance: transient growth of streamwise rollers

We attempt to provide an explanation of the weakening of VLSMs by orchestrating numerical experiments for the inspection of the generation of the very-large-scale velocity streaks excited by the streamwise rollers, which are the optimal output and input (Hwang & Cossu Reference Hwang and Cossu2010; McKeon Reference McKeon2017) or the optimal amplification and initial perturbations (Del Alamo & Jimenez Reference Del Alamo and Jimenez2006; Brandt Reference Brandt2014) obtained in linear stability analysis in wall-bounded turbulence (Alizard et al. Reference Alizard, Pirozzoli, Bernardini and Grasso2015; Chen et al. Reference Chen, Cheng, Fu and Gan2023). Such a task requires the accurate distribution of the mean velocity and an initial perturbation. The former is obtained via the postulation that (i) the mean velocities are consistent with those of the incompressible flows, independent of Mach numbers and wall temperatures, by applying the following transformation to the mean velocity and the wall-normal coordinate proposed by Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020):

(3.1a,b)\begin{equation} y^+_{VP} = \frac{1}{\delta_\nu} \int^y_0 \sqrt{\frac{\bar{\rho}}{\bar{\rho}_w} \frac{\bar{\mu}^3_w}{\bar{\mu}^3}} \,{\rm d} y,\quad u^+_{VP}= \frac{1}{u_\tau} \int^{\tilde u}_0 \left( \sqrt{\frac{\bar{\rho}}{\bar{\rho}_w} \frac{\bar{\mu}_w}{\bar{\mu}}}\right)\, {\rm d} \tilde u, \end{equation}

and that (ii) the mean velocity of the incompressible turbulent boundary layers can be constructed as (Subrahmanyam, Cantwell & Alonso Reference Subrahmanyam, Cantwell and Alonso2022)

(3.2)\begin{equation} \bar{u}^+(y^+) = \int^{y^+}_0 \left[ - \frac{1}{2\lambda (s)^2} + \frac{1}{2\lambda (s)^2} \left( 1+ 4 \lambda(s)^2 \left( 1 - \frac{s}{Re_\tau} \right)\right)^{1/2} \right]\, {\rm d}s, \end{equation}

in which the mixing-length function is expressed as

(3.3)\begin{equation} \lambda = \frac{k y^+ (1-\exp(-(y^+{/}a)^m))}{(1+(y^+{/}(b Re_\tau))^n)^{1/n}} \end{equation}

with $k=0.42$, $a=24.96$, $m=1.15$, $b=0.18$ and $n=2.17$. The total shear stress distribution is presumed to satisfy

(3.4a,b)\begin{equation} \tau(y)= { \bar{\mu}} \frac{\partial \bar{u}}{\partial y} + R_{uv}= \tau_w (1-\chi^{3/2}),\quad \chi = \frac{y/\delta}{(1+(y/\delta)^p)^{1/p}}, \end{equation}

with $p=17$ to ensure that $\tau$ gradually decays to zero as it reaches the free stream (Chen & She Reference Chen and She2016; Wang Reference Wang2022). With the aforementioned information and the generalized Reynolds analogy (Zhang et al. Reference Zhang, Bi, Hussain and She2014), the mean velocity can be readily obtained by setting the Mach number, wall temperature and the target friction Reynolds number. The constructed mean velocity and the Reynolds shear stress at $Re_\tau = 1000$, $M_\infty$ ranging from $0.3$ to $6.0$ and $T_w=T_r$, as shown in figure 4, are well-collapsed with those reported by Bernardini & Pirozzoli (Reference Bernardini and Pirozzoli2011), particularly that at $M_\infty =2.0$, suggesting the validity of the presently adopted methodology. Such a construction of the mean velocity and the Reynolds shear stress also shows considerable accuracy at lower Reynolds numbers, various Mach numbers and wall temperatures (Zhang, Duan & Choudhari Reference Zhang, Duan and Choudhari2018), which are omitted here for brevity.

Figure 4. Wall-normal distribution of constructed (a) $u^+_{VD}$ and (b) $-R^+_{uv}$. Symbols in panel (a): Bernardini & Pirozzoli (Reference Bernardini and Pirozzoli2011) at $M_\infty =2$ and $Re_\tau =1000$.

According to the study by Alizard et al. (Reference Alizard, Pirozzoli, Bernardini and Grasso2015) and Dawson & McKeon (Reference Dawson and McKeon2020), the optimal initial perturbations are weakly dependent on Mach numbers. Therefore, the initial perturbation in the form of streamwise vortices is generated using the following cross-stream velocity:

(3.5)\begin{equation} \left. \begin{aligned} v' & = A \exp(-(y/\delta)^2) \cos(2 {\rm \pi}z) \sin(2 {\rm \pi}\sqrt{y})/\sqrt{\bar{\rho}},\\ w' & = A \exp(-(y/\delta)^2) \sin(2 {\rm \pi}z) \sin(2 {\rm \pi}\sqrt{y})/\sqrt{\bar{\rho}}, \end{aligned} \right\} \end{equation}

with the initial amplitude $A=0.01 U_\infty$ and the mean density divided so as to ensure that the turbulent kinetic energies of the initial perturbation are the same for cases at different Mach numbers.

The numerical simulations of this process are performed in rectangular boxes with the sizes of $(10, 3, 1)\delta$ and the periodic conditions in the streamwise and spanwise directions to obtain the temporal evolution of the initial perturbation. The mean profiles are enforced to remain steady at each step. In figure 5(a), we present the turbulent kinetic energy amplification of the streamwise velocity $u'$ and cross-stream velocities $v'$ and $w'$, denoted by $G_{u}(t)$ and $G_{vw}(t)$, which are calculated by the ratio between the integrated energy of the corresponding velocity components in the computational domain and that of initial perturbations. For all the cases considered, the streamwise rollers (figure 5b, vectors) decay at the same rate, as can be inferred from their almost identical temporal evolution. These structures bring the low-speed fluids upwards and the high-speed fluids downwards, inducing the spanwise alternating structures resembling the VLSMs in moderate- and high-Reynolds-number flows. Intriguingly, the amplification of these large-scale streamwise velocity fluctuations $G_{u}(t)$ decreases with the Mach number, despite that the kinetic energy of initial perturbations is determined to be identical, and that the optimal growth rates are reached at the same time instant. Therefore, it can be concluded that the abatement of the streamwise velocity fluctuations should be attributed to the inherent properties of the flow dynamics that the perturbations of the same energy are only capable of triggering the VLSM structures with lower energy at higher Mach numbers, in accordance with the statistical results in § 2. The conclusion is the same by setting the initial perturbation with the same magnitude of velocity (without weighted by density in (3.5)) and momentum (divided by $\bar {\rho }$ in (3.5)). The results obtained in this way are qualitatively the same as those obtained in compressible turbulent channels by Chen et al. (Reference Chen, Cheng, Fu and Gan2023) under the framework of linearized governing equations with harmonic and stochastic forcing.

Figure 5. (a) Energy amplification of the initial perturbation in the form of streamwise vortices: solid, $G_{vw}(t)$; dotted, $G_{u}(t)$; line legends refer to figure 4. (b) Initial $v'$ and $w'$ perturbation (vector) and the $u'$ at $t=9 \delta /U_\infty$ (flooded).

4. Turbulent kinetic energy production

Inspecting the turbulent production term that extracts the kinetic energy from the mean flow to turbulent fluctuations is another possibly effective perspective for explaining the variation of turbulent kinetic energy at different flow parameters. One example would be that the increasing turbulent production in the outer region corresponds to the enhanced turbulent kinetic energy and prominent VLSMs centred at $y \approx 0.2\delta$ and that a long extending plateau in the logarithmic layer is manifested at sufficiently high Reynolds numbers (Marusic, Mathis & Hutchins Reference Marusic, Mathis and Hutchins2010; Smits et al. Reference Smits, McKeon and Marusic2011). According to its expression $P_K = R_{uv} {\partial \bar {u}}/{\partial y}$, the turbulent production term can be constructed once the mean velocity and the Reynolds shear stress are obtained from (3.1a,b)–(3.4a,b).

In figure 6(a), we present the distribution of $P_K$ normalized by viscous scales and pre-multiplied by $y^+$ for a better representation of the integration against the logarithmic abscissa. The peaks in the outer region located at $y^+ \approx 400$ decrease monotonically to be manifested as a plateau at $M_\infty =4.0$ and as a shoulder at $M_\infty > 5.0$. The conclusions remain valid when inspected by the LES results in § 2. This confirms our previous elucidations and postulations that the turbulent intensities are weakened by the higher Mach numbers. Such a variation can also be inferred by reformulating the production term as

(4.1)\begin{equation} P^+_K = R^+_{uv} \frac{\partial \bar{u}^+}{\partial y^+} = R^+_{uv} \frac{\partial \bar{u}^+_{VD}}{\partial y^+} \sqrt{\frac{\bar{\rho}_w}{\bar{\rho}}} \approx R^+_{uv} \frac{\partial \bar{u}^+_{VD}}{\partial y^+} \sqrt{\frac{\bar{T}}{ T_w}}, \end{equation}

where the mean pressure across the boundary layer is presumed to be constant. Since we have shown in figure 4 that the van Driest transformed mean velocity $\bar {u}^+_{VD}$ and the Reynolds shear stress $R^+_{uv}$ are weakly dependent on the Mach number, the first two factors in the last expression of (4.1) should be regarded as invariants, leaving mean temperature variation as the only source of disparity. For high-speed turbulent boundary layers over adiabatic walls ($T_w = T_r$), it is well established that the ratio $\bar {T}/T_w$ decreases with the increment of Mach numbers at a certain off-wall location within the boundary layer. Henceforth, it is straightforward that the production term $P^+_K$, the rate of energy transfer from the mean to turbulent kinetic energy, is curtailed. Naturally, the profiles of the $P^+_K$ at various Mach numbers can be collapsed under the transformation proposed by Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020), but this means nothing, for the integration in the physical space ($y/\delta$) is identical, irrelevant of the transformation adopted.

Figure 6. Wall-normal distribution of the pre-multiplied constructed turbulent production $y^+ P^+_k$ at different (a) Mach numbers ($Re_\tau =1000$, $T_w=T_r$), (b) Reynolds numbers ($T_w=T_r$, dashed, $M_\infty =0.5$; solid, $M_\infty =6.0$) and (c) wall temperatures ($Re_\tau = 2000$, $M_\infty =6.0$), and (d) $y^+_V P^*_{K}$ normalized by semi-local scales, parameters are the same as those in panel (c). Symbols in panel (a): Bernardini & Pirozzoli (Reference Bernardini and Pirozzoli2011) at $M_\infty =2$ and $Re_\tau =1000$.

The method of construction of the turbulent production $P^+_K$ allows the exploration of turbulent characteristics without performing formidable and extravagant high-precision numerical simulations, especially when the Reynolds numbers are extremely high. Taking its advantage, we further discuss the influences of the other two parameters, i.e. the friction Reynolds number $Re_\tau$ and the wall temperature $T_w$.

In figure 6(b), we present the pre-multiplied turbulent production $y^+ P^+_k$ at the Mach numbers of $M_\infty =0.5$ and $6.0$, and the wall temperature of $T_w=T_r$ (adiabatic), with the friction Reynolds number ranging from $Re_\tau =500$ to $20\,000$. At the lower Mach number $M_\infty =0.5$, the variation of $y^+ P^+_k$ resembles that of the incompressible turbulent boundary layers in that the inner and outer peaks are gradually separated by a gradually elongated plateau in the logarithmic region. When it comes to high-Mach-number flows, however, the peaks and the plateau exist no more, merely manifesting longer ranges of logarithmic decay and small shoulder before they decrease abruptly to insignificant values. This suggests that the lower intensity of large- and very-large-scale turbulent structures and the turbulent intensities in the outer region are universal phenomena at all Reynolds numbers instead of a sporadic occurrence at certain combinations of flow parameters.

We also exploit the influences of the wall temperature by inspecting the production term $P^+_K$ at $M_\infty =6$ and $Re_\tau =2000$ with the wall temperature being $T_w=(1.0, 0.8, 0.6, 0.4, 0.2) T_r$. The results are shown in figure 6(c,d), normalized by viscous and semi-local scales, respectively. Expectedly, the turbulent production $P^+_K$ is well collapsed under the semi-local scaling (figure 6d) in the near-wall region and shows a long extending plateau in the log layer, resembling those of the incompressible flows. Nevertheless, such rescaled wall-normal coordinates are non-physical, as pointed out previously, for they are not uniformly transformed from the physical space, leading to the misconception of the wall temperature independence. When plotted against the viscous scales, as demonstrated in figure 6(c), the $P^+_K$ is augmented by the decreasing wall temperature, indicating that the cooling walls tend to enhance the large-scale motions in the outer region, at any rate in the perspective of kinetic energy transfer. However, the intervals between the two peaks (or the shoulders) are shortened, consistent with the DNS results given by Huang et al. (Reference Huang, Duan and Choudhari2022), suggesting the smaller scale separation between the inner and outer motions (Smits et al. Reference Smits, McKeon and Marusic2011). This has been verified to be consistent with the DNS results reported by Zhang et al. (Reference Zhang, Duan and Choudhari2018) at $M_\infty \approx 6$ and $Re_\tau =450$. The same conclusions can be drawn if transient growth of the large-scale streamwise rollers is used to investigate this issue, which is demonstrated in the Appendix.

To summarize, over adiabatic walls, the turbulent production term that represents the kinetic energy transfer from the mean to turbulent flows decreases monotonically with the increasing Mach number at all Reynolds numbers in the logarithmic and outer layers. At the same $Re_\tau$ and $M_\infty$, the turbulent productions above the buffer region are gradually increased by the cooling walls. Although they only serve as qualitative indicators of the intensification or diminishment of VLSMs, these elucidations are consistent with the existing DNS or LES results.