Principle of fundamental resonance in hypersonic boundary layers: an asymptotic viewpoint

Abstract

The fundamental resonance (FR) in the nonlinear phase of the boundary-layer transition to turbulence appears when a dominant planar instability mode reaches a finite amplitude and the low-amplitude oblique travelling modes with the same frequency as the dominant mode, together with the stationary streak modes, undergo the strongest amplification among all the Fourier components. This regime may be the most efficient means to trigger the natural transition in hypersonic boundary layers. In this paper, we aim to reveal the intrinsic mechanism of the FR in the weakly nonlinear framework based on the large-Reynolds-number asymptotic technique. It is found that the FR is, in principle, a triad resonance among a dominant planar fundamental mode, a streak mode and an oblique mode. In the major part of the boundary layer, the nonlinear interaction of the fundamental mode and the streak mode seeds the growth of the oblique mode, whereas the interaction of the oblique mode and the fundamental mode drives the roll components (transverse and lateral velocity) of the streak mode, which leads to a stronger amplification of the streamwise component of the streak mode due to the lift-up mechanism. This asymptotic analysis clearly shows that the dimensionless growth rates of the streak and oblique modes are the same order of magnitude as the dimensionless amplitude of the fundamental mode $(\bar {\epsilon }_{10})$, and the amplitude of the streak mode is $O(\bar {\epsilon }_{10}^{-1})$ greater than that of the oblique mode. The main-layer solution of the streamwise velocity, spanwise velocity and temperature of both the streak and the oblique modes become singular as the wall is approached, and so a viscous wall layer appears underneath. The wall layer produces an outflux velocity to the main-layer solution, inclusion of which leads to an improved asymptotic theory whose accuracy is confirmed by comparing with the calculations of the nonlinear parabolised stability equations (NPSEs) at moderate Reynolds numbers and the secondary instability analysis (SIA) at sufficiently high Reynolds numbers.

1. Introduction

In the development of hypersonic vehicles, to calculate accurately the drag and heat flux is a task of the first priority, which requires an accurate prediction of the laminar–turbulent transition. For high-altitude flight conditions, transition is often triggered by a natural route, for which four phases, including the receptivity, linear instability, nonlinear resonance and turbulence, appear in sequence. For supersonic or hypersonic boundary layers, there may exist more than one discrete instability mode, which are referred to as the Mack first, second, $\ldots$ modes, according to the ascending order of their frequencies (Mack 1987). It was revealed by the asymptotic analysis that only the Mack first mode with $\varTheta > \tan ^{-1} \sqrt {M^2-1}$ (where $\varTheta$ and $M$ denote the wave angle and Mach number, respectively) has the viscous nature (Smith 1989; Liu, Dong & Wu 2020), while the quasi-two-dimensional Mack first and all the higher-order modes are inviscid (Cowley & Hall 1990; Smith & Brown 1990; Dong, Liu & Wu 2020; Zhao, He & Dong 2023). Usually, the Mack second mode appears when the Mach number is approximately over 4, and its growth rate peaks when it is planar (two-dimensional). In contrast, the Mack first mode appears in all supersonic boundary layers, which is more unstable when it is oblique (three-dimensional). The linear evolution of these modes was confirmed by quite a few numerical works, such as Fedorov (2011) and Zhong & Wang (2012). When the unstable modes are accumulated to finite amplitudes, the nonlinear interaction among different Fourier components becomes the leading-order impact, showing three major nonlinear regimes, including the oblique-mode breakdown (Thumm 1991; Fasel, Thumm & Bestek 1993; Chang & Malik 1994; Leib & Lee 1995; Mayer, Von Terzi & Fasel 2011a; Mayer, Wernz & Fasel 2011b), the subharmonic resonance (Saric 1984; Herbert 1988) and the fundamental resonance (Sivasubramanian & Fasel 2015; Hader & Fasel2019).

The oblique-mode breakdown appears when the dominant perturbations in the early nonlinear phase are a pair of three-dimensional (3-D) travelling waves with the same frequency but opposite spanwsie wavenumbers, as sketched in figure 1(a). Such a regime was pioneered by Thumm (1991) and Fasel et al. (1993) using the direct numerical simulation (DNS) approach and subsequently studied by Chang & Malik (1994) using the nonlinear parabolised stability equations (NPSEs) approach. They all found that the growth rates of the second-order harmonics and the stationary streak mode are equal to twice those of the oblique modes, because they are driven by the self or mutual interaction of the oblique modes. Particularly, the streak mode shows a much greater amplitude than the second-order harmonics. Such a phenomenon was recently explained by Song, Dong & Zhao (2022) using the weakly nonlinear analysis based on the large-Reynolds-number asymptotic technique. Considering the growth rate of the travelling Mack mode to be much smaller than its wavenumber, the transverse and lateral perturbation velocities of the streak mode, showing a roll structure, are primarily amplified due to the mutual interaction of the oblique modes, but its streamwise velocity perturbation undergoes a further amplification due to the lift-up mechanism induced by the roll structure. For a low-Mach-number supersonic boundary layer, the most unstable linear perturbation is usually the oblique first mode, which ensures the dominant perturbation in the early nonlinear phase to be three-dimensional, indicating that the oblique-mode breakdown regime is likely to be triggered in this configuration.

Figure 1. Characteristic parameters in the frequency-spanwise-wavenumber space: (a) oblique-mode breakdown regime; (b) subharmonic resonance regime; (c) fundamental resonance regime. Here, $\omega _0$ and $\beta _0$ denote the fundamental frequency and spanwise wavenumber, respectively. This sketch is a modified version of figure 10 of Hader & Fasel (2019).

The subharmonic resonance (SR) appears when the dominant perturbation in the early nonlinear phase is planar, or two-dimensional (2-D), and the frequency of the most amplified 3-D perturbations is half of that of the 2-D mode; see the schematic in figure 1(b). Such a regime usually appears in a subsonic or an incompressible boundary layer, as observed numerically (Herbert 1988) and experimentally (Saric 1984), and the rapid amplification of the 3-D modes is attributed to the secondary instability (SI) based on the Floquet theory (Herbert 1988). For supersonic boundary layers, Kosinov, Maslov & Shevelkov (1990), Kosinov et al. (1994) and Kosinov & Tumin (1996) reported a generalised subharmonic regime, for which the dominant perturbation is a 3-D Mack first mode, and the two most promoted SI modes are subharmonic in frequency. This scenario was later confirmed by the numerical simulations of Mayer et al.(2011b).

For a hypersonic boundary layer, the 2-D Mack second mode is the most linearly amplified perturbation, which could be the dominant perturbation in the early nonlinear phase. Using the DNS approach, the most amplified 3-D modes are found to be those with the same frequency as the 2-D second mode (Sivasubramanian & Fasel 2015; Hader & Fasel 2019), which was also confirmed by the secondary instability analysis (SIA) of Chen, Zhu & Lee (2017). This scenario is referred to as the fundamental resonance (FR) and a schematic for the FR in the spectrum space is shown in figure 1(c). Based on the critical-layer theory (Wu 2004; Zhang & Wu 2022), Wu, Luo & Yu (2016) deduced the evolution equations for the oblique modes and claimed that the 2-D mode acts as a catalyst to promote the growth of the oblique modes. The 2-D dominant mode and the small-amplitude oblique modes are found to be phase-locked. Actually, the SI modes include both the 3-D travelling waves and the stationary streak mode, and the amplitude of the latter was found to be much greater than those of the former (Brad et al. 2009; Chou et al. 2011; Sivasubramanian & Fasel 2015; Chynoweth et al. 2019; Hader & Fasel 2019). However, the latter phenomenon so far is not well explained from the dynamicviewpoint.

For convenience of illustration, each Fourier component with a frequency $m\omega _0$ and a spanwise wavenumber $n\beta _0$ is denoted by $(m,n)$, where $\omega _0$ and $\beta _0$ are the fundamental frequency and spanwise wavenumber, respectively. It is seen that in both the SR and FR regimes, the dominant perturbation in the early nonlinear phase is 2-D. If we choose the frequency of the 2-D mode to be the fundamental frequency $\omega _0$, then the 2-D fundamental mode is denoted by $(1,0)$. For the SR, the most unstable 3-D travelling waves are components $(1/2,n)$, where $n$ is an integer to represent the spanwise wavenumber. The components $(1,0)$, $(1/2,n)$ and $(1/2,-n)$ with a non-zero $n$ form a triad resonance system, for which the mutual interaction of any two components seeds for the growth of the third one. When the dominant 2-D mode $(1,0)$ reaches a nonlinear saturation phase, the oblique modes $(1/2,n)$ and $(1/2,-n)$ could amplify with a greater rate, interpreted as an SI regime. As the oblique pair $(1/2,\pm n)$ reach finite amplitudes, their interaction could also drive a streak mode $(0,2n)$, similar to the oblique-mode breakdown regime. However, for the FR, the most unstable 3-D travelling waves are $(1,n)$ with $n$ being an integer, and the Fourier components $(1,0)$, $(1, n)$ and $(1,-n)$ with a non-zero $n$ do not form a triad resonance system. Although SI analyses have confirmed the rapid amplification of $(1,\pm n)$ in the nonlinear phase (Sivasubramanian & Fasel 2015; Chen et al. 2017; Hader & Fasel 2019), the energy transfer among different Fourier components due to their nonlinear interaction is far beyond obvious. Actually, the SI modes include a set of oblique modes with the same frequency $(1,n)$ and stationary streak modes $(0,n)$, and the Fourier components $(1,0)$, $(1,n)$ and $(0,n)$ could form a triad resonance system. Inclusion of the streak mode in the triad resonance determines the key role of the streak mode in the SI process of the FR, which is in contrast to the SR. Unfortunately, such a dynamic mechanism has not been formulated theoretically, especially in the hypersonic boundary layers, which is the main task of the present work.

The rest part of this paper is structured as follows. In § 2, we introduce the physical model and the numerical treatment (NPSE approach) for the FR, and the NPSE calculations showing the crucial role of the streak mode are demonstrated in § 3. In § 4, we develop an asymptotic theory to describe the dynamic mechanism of the FR, whose accuracy is confirmed by the NPSE calculations for moderate Reynolds numbers in § 5 and by SI analysis for sufficiently high Reynolds numbers in § 6. The concluding remarks are present in § 7.

2. Physical model and governing equations

2.1. Physical model

The physical model to be studied is a flat plate inserted into a perfect-gas hypersonic stream with zero angle of attack, as sketched in figure 2. The plate is assumed to be infinitely thin such that a rather weak oblique shock forms from its leading edge, and the flow quantities behind the shock are rather close to those of the oncoming stream. A viscous boundary layer forms in the close neighbourhood of the wall.

Figure 2. Sketch of the physical model.

We study the evolution of a set of Mack instability modes in a selected computational domain, whose inlet boundary is located at a distance $L^*$ downstream of the leading edge, where the Mack modes are already unstable. The inflow perturbations consist of a dominant 2-D Mack mode and a pair of small 3-D Mack modes with the same frequency. The problem is described by the Cartesian coordinate $\boldsymbol {x}^*=(x^*,y^*,z^*)$ with its origin $O$ at the inlet of the computational domain. The reference length is selected to be the boundary-layer characteristic thickness at the origin $\delta ^*=\sqrt {\nu ^{*}_\infty L^{*}/U^{*}_\infty }$, where $U_\infty ^*$ and $\nu ^{*}_\infty$ denote the velocity and the kinematic viscosity of the oncoming flow. The dimensionless coordinate and time are expressed as

(2.1a,b)\begin{equation} \boldsymbol{x}=(x,y,z)=\boldsymbol{x}^*/\delta^*,\quad t=t^*U^{*}_\infty/\delta^*. \end{equation}

In what follows, the superscript $*$ denotes the dimensional quantities. The velocity field $\boldsymbol {u}=(u,v,w)$, density $\rho$, temperature $T$, pressure $p$ are normalised by $U^{*}_\infty$, $\rho ^{*}_\infty$, $T^{*}_\infty$, $\rho ^{*}_\infty U^{{*}2}_\infty$, respectively, where the subscript $\infty$ denotes the quantities at the oncoming stream. For unsteady perturbations, the frequency $\omega$, streamwise wavenumber $\alpha$ and spanwise wavenumber $\beta$ are normalised as

(2.2a–c)\begin{equation} \omega=\omega^* \delta^*/U^{*}_\infty,\quad \alpha=\alpha^* \delta^*,\quad\beta=\beta^*\delta^*. \end{equation}

The flow field is governed by two dimensionless parameters, the Reynolds number $R=U^{*}_\infty \delta ^*/\nu ^{*}_\infty$ and the Mach number $M=U_\infty ^*/a_\infty ^*$, with $a_\infty ^*$ denoting the sound speed of the oncoming stream.

2.2. Governing equations

The dimensionless governing equations are

(2.3a)$$\begin{gather} \frac{{\partial \rho }}{{\partial t}} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( {\rho {\boldsymbol{u}}} \right) = 0, \end{gather}$$

(2.3b)$$\begin{gather}\rho {\frac{{{\rm D}{\boldsymbol{u}}}}{{{\rm D}t}} } ={-} \boldsymbol{\nabla} {p} + \frac{1}{{{{R}}}}\left[ {2\boldsymbol{\nabla} \boldsymbol{\cdot} \left( {\mu {\boldsymbol{S}}} \right) - \frac{2}{3}\boldsymbol{\nabla} \left( {\mu \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}}} \right)} \right], \end{gather}$$

(2.3c)$$\begin{gather}\frac{1}{\gamma }\rho \frac{{{\rm D}T}}{{{\rm D}t}} - \frac{{\gamma - 1}}{\gamma }T\frac{{{\rm D}\rho }}{{{\rm D}t}} = \frac{1}{{PrR}}\boldsymbol{\nabla} \boldsymbol{\cdot} \left( {\mu \boldsymbol{\nabla} T} \right) + \frac{{\left( {\gamma - 1} \right){M^2}}}{{R}}\left[ {2\mu {\boldsymbol{S}}:{\boldsymbol{S}} - \frac{2}{3}\mu {{\left( {\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}}} \right)}^2}} \right], \end{gather}$$

(2.3d)$$\begin{gather}{p} = \frac{{\rho }{T}}{\gamma {M}^2}, \end{gather}$$

where ${\boldsymbol {S}} = {{[ {\boldsymbol {\nabla } {\boldsymbol {u}} + {{( {\boldsymbol {\nabla } {\boldsymbol {u}}} )}^{\rm T}}} ]}/ 2}$ is the strain-rate tensor, $\mu$ is the dynamic viscosity satisfying the Sutherland law ($\mu =(1+ { C})T^{{3}/{2}}/(T+ { C})$ with ${C}=110.4K/{T_\infty ^{*}}$), $Pr$ is the Prandtl number, $\gamma$ is the ratio of the specific heats and ${{\rm D}}/{{\rm D}t} = {\partial }/{\partial t} + \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla }$ denotes the material derivative. In this paper, we take $Pr=0.72$ and $\gamma =1.4$.

The no-slip, non-penetration and isothermal boundary conditions are applied at the wall,

(2.4)\begin{equation} (u,v,w,T)=(0,0,0,T_w), \quad \text{at} \ y=0, \end{equation}

where $T_w$ is the dimensionless wall temperature. In the far field, all the perturbations damp exponentially (the radiating mode as that of Chuvakhov & Fedorov (2016) is not considered here) and the upper boundary conditions read

(2.5)\begin{equation} (\rho, u,v,w,T)\rightarrow(1,1,0,0,1), \quad \text{as} \ y \to \infty. \end{equation}

2.3. Perturbation field

The instantaneous flow field $\phi \equiv (\rho,u,v,w, T)$ can be decomposed into a steady base flow $\varPhi _B$ and an unsteady perturbation $\tilde {\phi }$,

(2.6)\begin{equation} \phi = \varPhi_B (x,y) + \tilde{\phi} (x,y,z,t), \end{equation}

where $\varPhi _B=(1/T_B,U_B,V_B,0,T_B)$ is the compressible Blasius solution. Substituting (2.6) into (2.3), and subtracting the base flow out, we obtain the nonlinear equations governing the perturbations,

(2.7)$$\begin{gather} {\boldsymbol { G }}\frac{{\partial \tilde{\phi} }}{{\partial t}} + {\boldsymbol{A}}\frac{{\partial \tilde{\phi} }}{{\partial x}} + {\boldsymbol{B}}\frac{{\partial \tilde{\phi} }}{{\partial y}} + {\boldsymbol{C}}\frac{{\partial \tilde{\phi} }}{{\partial z }} + {\boldsymbol{D}}\tilde{\phi} + {{\boldsymbol{V}}_{xx}}\frac{{{\partial ^2}\tilde{\phi} }}{{\partial {x^2}}} + {{\boldsymbol{V}}_{yy}}\frac{{{\partial ^2}\tilde{\phi} }}{{\partial {y^2}}} \nonumber\\ +\, {{\boldsymbol{V}}_{z z }}\frac{{{\partial ^2}\tilde{\phi} }}{{\partial {z ^2}}} + {{\boldsymbol{V}}_{xy}}\frac{{{\partial ^2}\tilde{\phi} }}{{\partial x\partial y}} + {{\boldsymbol{V}}_{xz }}\frac{{{\partial ^2}\tilde{\phi} }}{{\partial x\partial z }} + {{\boldsymbol{V}}_{yz }}\frac{{{\partial ^2}\tilde{\phi} }}{{\partial y\partial z }} = \boldsymbol{F}, \end{gather}$$

where the coefficient matrices $\boldsymbol {G}$, $\boldsymbol {A}$, $\boldsymbol {B}$, $\boldsymbol {C}$, $\boldsymbol {D}$, ${\boldsymbol {V}}_{xx}$, ${\boldsymbol {V}}_{yy}$, ${\boldsymbol {V}}_{z z}$, ${\boldsymbol {V}}_{xy}$, ${\boldsymbol {V}}_{yz}$ and ${\boldsymbol {V}}_{x z}$ and the nonlinear forcing $\boldsymbol {F}$ can be found in Appendix A. The pressure perturbation $\tilde {p}$ has been eliminated by the equation of the state (2.3d).

In this paper, we particularly focus on the FR regime, therefore, the introduced perturbations $\tilde {\phi }$ at the inlet of the computational domain should include a 2-D Mack second mode with a frequency $\omega _0$ and a pair of 3-D Mack second modes with the same frequency $\omega _0$ but opposite spanwise wavenumbers $\pm \beta _0$, where $\omega _0$ and $\beta _0$ are referred to as the fundamental frequency and fundamental spanwise wavenumber, respectively. For convenience, we use $(m,n)$ to denote a perturbation with a frequency $m\omega _0$ and a spanwise wavenumber $n\beta _0$. Thus, the three introduced perturbations are denoted by $(1,0)$, $(1,1)$ and $(1,-1)$, respectively. The inflow perturbation can be expressed as

(2.8)$$\begin{align} \tilde{\phi}(0,y,z,t) &=\epsilon_{10}\hat \phi_{10}(\kern0.7pt y)\exp(-{\rm i}\omega_0 t)+\epsilon_{11}\hat \phi_{11}(\kern0.7pt y)\exp[{\rm i}(\beta_0z-\omega_0 t)] \nonumber\\ &\quad +\epsilon_{1-1}\hat \phi_{1-1}(\kern0.7pt y)\exp[{\rm i}(-\beta_0z-\omega_0t)]+ \text{c.c.}, \end{align}$$

where $\epsilon _{mn}$ measures the initial amplitude of the introduced perturbation $(m,n)$, $\hat {\phi }_{mn}$ denotes the perturbation profile for the $(m,n)$ component, $\text {c.c.}$ denotes the complex conjugation and $\text {i} \equiv \sqrt {-1}$. For a hypersonic boundary layer, the 2-D second mode is usually more unstable, and the amplitude of mode $(1,0)$ should be much greater than modes $(1,\pm 1)$ due to the historical accumulative effect. Thus, the amplitude of the introduced 2-D mode $\epsilon _{10}$ is taken to be much greater than those of the 3-D modes $\epsilon _{11}$ and $\epsilon _{1-1}$. Also, the linear growth rates of the two oblique modes are the same, and so we let $\epsilon _{11}=\epsilon _{1-1}$. Note that in reality, the spanwise wavenumbers of the oblique modes may be not exactly opposite and their amplitudes may be different, but our selection is still a good demonstration to reveal their resonance mechanism.

2.3.1. Linear stability theory (LST)

The perturbation profile $\hat \phi$ for each Fourier mode of the inflow perturbation is obtained by the linear stability theory (LST). Introducing the parallel-flow assumption, the perturbation with a frequency $\omega$, a streamwise wavenumber $\alpha$ and a spanwise wavenumber $\beta$ is expressed in terms of a travelling-wave form,

(2.9)\begin{equation} \tilde{\phi} = \epsilon_L \hat{\phi}(\kern0.7pt y) \exp [{\rm i} (\alpha x + \beta z -\omega t) ] + \text{c.c.}, \end{equation}

where $\epsilon _L \ll 1$ measures its amplitude. Substituting (2.9) into the system (2.7) with $O(\epsilon _L^{2})$ terms neglected, we obtain the compressible Orr–Sommerfeld (O-S) equations,

(2.10)\begin{equation} {\boldsymbol{\tilde{B}}}\frac{{\partial \hat \phi }}{{\partial y}} + {{\boldsymbol{V}}_{yy}}\frac{{{\partial ^2}\hat \phi }}{{\partial {y^2}}} + {\boldsymbol{\tilde{D}}}\hat \phi = 0, \end{equation}

where

(2.11)\begin{equation} \left. \begin{gathered} {\boldsymbol{\tilde{B}}} = {\boldsymbol{B}} + {\rm{i}}\alpha {{\boldsymbol{V}}_{xy}} + {\rm{i}}\beta{{\boldsymbol{V}}_{yz }},\\ {\boldsymbol{\tilde{D}}} ={-} {\rm{i}}\omega {\boldsymbol{G}} + {\rm{i}}\alpha {\boldsymbol{A}} + {\rm{i}}\beta{\boldsymbol{C}} + {\boldsymbol{D}} - {\alpha ^2}{{\boldsymbol{V}}_{xx}} - {\beta^2}{{\boldsymbol{V}}_{zz }} - \alpha \beta{{\boldsymbol{V}}_{xz }}. \end{gathered} \right\} \end{equation}

Introducing the homogeneous boundary conditions,

(2.12a,b)\begin{equation} [\hat{u},\hat{v},\hat{w},\hat{T}]=0, \quad \text{at} \ y=0; \quad \hat{\phi} \to 0, \quad \text{as} \ y\to \infty, \end{equation}

we arrive at an eigenvalue problem. For the spatial mode, $\omega$ and $\beta$ are given to be real, and the eigenvalue $\alpha = \alpha _r + \textrm {i} \alpha _i$ is complex with the opposite of its imaginary part representing the growth rate. Usually, the imaginary part is much smaller than the real part in the boundary-layer flow, i.e. $|\alpha _i| \ll |\alpha _r|$. The numerical details to solve the eigenvalue system (2.10) with (2.12a,b) can be found in our previous papers (Dong et al. 2020; Song, Zhao & Huang 2020; Dong & Zhao 2021; Li & Dong 2021).

2.3.2. Nonlinear parabolised stability equations (NPSEs)

The NPSE approach (Bertolotti, Herbert & Spalart 1992; Chang & Malik 1994) is considered as a more accurate means because it allows the slow streamwise variation of the perturbation profiles and takes into account the non-parallelism of the base flow. The only approximation is that the $\partial _{xx}$ terms are neglected to reduce the elliptic system to a parabolised system, which is quite reasonable for a boundary layer with a smooth wall. Expressing $\tilde {\phi }$ and $\boldsymbol {F}$ in terms of the Fourier series with respect to $z$ and $t$, we obtain

(2.13)\begin{equation} \left. \begin{gathered} \tilde{\phi}(x,y,z,t) = \sum_{m ={-}M_e}^{{M_e}} {\sum_{n ={-} {N_e}}^{{N_e}} {\mathop{\phi}\limits^{{\smile}}}_{mn}(x,y)\exp [\text{i} \left(n {\beta_0} z - m {\omega_0} t\right) ] } ,\\ {\boldsymbol{F}} (x,y,z,t) = \sum_{m ={-}M_e}^{{M_e}} {\sum_{n ={-} {N_e}}^{{N_e}} {{{{\boldsymbol{\tilde{F}}}}_{mn} (x,y) }{ \exp [ {\text{i}\left( {n{\beta_0}z - m{\omega _0}t} \right)} ] }} } , \end{gathered} \right\} \end{equation}

where $M_e$ and $N_e$ denote the orders of the Fourier-series truncation. In this paper, we choose $M_e =5$ and $N_e =5$, which has been confirmed to be sufficient via resolution tests. Considering that the perturbations are propagating with two length scales, a fast one with an oscillatory manner and a slow one related to the non-parallelism, we express the perturbation profile $\mathop {\phi }\limits ^{\smile }$ in terms of a Wentzel–Kramers–Brillouin (WKB) form,

(2.14)\begin{equation} {{\mathop{\phi}\limits^{{\smile}}}_{mn}} (x,y) = {{\check \phi }_{mn}}\left( {x,y} \right){\exp\left({{\rm i} \int_{{x_0}}^x {{\alpha _{mn}} (\bar x) \,{\rm d}\kern0.06em \bar x} }\right)}, \end{equation}

where each Fourier component is denoted by $(m,n)$, $\omega _{0}$ and $n_{0}$ are the fundamental frequency and spanwise wavenumber, respectively, and $\alpha _{mn}$ represents the complex streamwise wavenumber of $(m,n)$. The shape function $\check {\phi }_{mn}$ varies slowly with $x$. The integral in (2.14) starts from a reference streamwise position $x_0$, which is selected as the inlet of the computational domain for the numerical calculations in this paper, namely, $x_0 \equiv 0$.

Neglecting the $\partial _{xx} \check {\phi }_{mn}$ terms, (2.7) is reduced to

(2.15)\begin{equation} {{{\boldsymbol{\tilde{A}}}}_{mn}}\frac{{\partial {{\check \phi }_{mn}}}}{{\partial x}} + {{{\boldsymbol{\tilde{B}}}}_{mn}}\frac{{\partial {{\check \phi }_{mn}}}}{{\partial y}} + {{\boldsymbol{V}}_{yy}}\frac{{{\partial ^2}{{\check \phi }_{mn}}}}{{\partial {y^2}}} + {{{\boldsymbol{\tilde{D}}}}_{mn}}{{\check \phi }_{mn}} = {{{{{\boldsymbol{\check F}}}_{mn}}}}, \end{equation}

where the matrices ${{{\boldsymbol {\tilde {A}}}}_{mn}}$, ${{{\boldsymbol {\tilde {B}}}}_{mn}}$ and ${{{\boldsymbol {\tilde {D}}}}_{mn}}$ are given by

(2.16)\begin{equation} \left. \begin{aligned} {{{\boldsymbol{\tilde{A}}}}_{mn}} & = {\boldsymbol{A}} + 2{\rm{i}}{\alpha _{mn}}{{\boldsymbol{V}}_{xx}} + {\rm{i}}n{\beta_0}{{\boldsymbol{V}}_{xz }},\\ {{{\boldsymbol{\tilde{B}}}}_{mn}} & = {\boldsymbol{B}} + {\rm{i}}{\alpha _{mn}}{{\boldsymbol{V}}_{xy}} + {\rm{i}}n{\beta_0}{{\boldsymbol{V}}_{yz }},\\ {{{\boldsymbol{\tilde{D}}}}_{mn}} & ={-} {\rm{i}}m{\omega _0} {\boldsymbol {G}} + {\rm{i}}{\alpha _{mn}}{\boldsymbol{A}} + {\rm{i}}n{\beta_0}{\boldsymbol{C}} + {\boldsymbol{D}}- {n^2}\beta_0^2{{\boldsymbol{V}}_{z z }}\\ & \quad - \left( {\alpha _{mn}^2 - {\rm{i}}\frac{{{\rm d}{\alpha _{mn}}}}{{{\rm d}\kern0.06em x}}} \right){{\boldsymbol{V}}_{xx}} - n{\alpha _{mn}}{\beta_0}{{\boldsymbol{V}}_{xz }},\\ {{{{{\boldsymbol{\check F}}}_{mn}}}} & = {{{{{\boldsymbol{\tilde{F}}}}_{mn}}}}{{{\exp({-\text{i}\int_{{x_0}}^x {{\alpha _{mn}} (\bar x)\,{\rm d}\kern0.06em \bar x} })}}}. \end{aligned} \right\} \end{equation}

The inflow perturbations are given by (2.8), and the lower and upper boundary conditions are

(2.17)\begin{equation} \left. \begin{gathered} ({{\check u}_{mn}},{{\check v}_{mn}},{{\check w}_{mn}},{{\check T}_{mn}}) = (0,0,0,0), \quad \text{at} \ y = 0 ,\\ ({{\check \rho}_{mn}},{{\check u}_{mn}},{{\check v}_{mn}},{{\check w}_{mn}},{{\check T}_{mn}}) \to (0,0,0,0,0), \quad \text{as} \ y \to \infty. \end{gathered} \right\} \end{equation}

To solve for the complex streamwise wavenumber $\alpha _{mn}$ and the profiles $\check {\phi }$, an iterative procedure is employed, which can be found from Zhao et al. (2016), and our code validation is provided in the appendix of Song et al. (2022).

Additionally, if $\check {\boldsymbol {F}}_{mn}$ is set to be zero, then (2.15) is recast to the linear PSE (LPSE), which can be used to track the evolution of each linear mode individually. To be distinguished, the PSE approach with $\check {\boldsymbol {F}}_{mn}$ being retained is referred to as the nonlinear PSE (NPSE) in this paper.

2.3.3. Secondary instability analysis (SIA) for a wavy base flow

When the amplitude of the fundamental perturbation $(1,0)$ has reached a finite level, the rapid growth of the infinitesimal perturbations can be explained by the SI based on a wavy profile driven by a 2-D quasi-saturated travelling mode.

Since the growth rate of the fundamental mode is usually much smaller than that of the SI mode, as confirmed by many numerical studies such as those of Sivasubramanian & Fasel (2015), Chen et al. (2017) and Hader & Fasel (2019), we take $\alpha _{10}$ to leading order to be real and introduce $\breve \alpha \equiv \mathrm {Re} (\alpha _{10})$. Thus, the base flow for the SIA is a superposition of the steady base flow and a series of quasi-saturated travelling waves. In a moving frame, the base flow $\breve {\varPhi }_{B} \equiv [\breve \rho _B,\breve U_B,0,0,\breve T_B]$ is expressed as

(2.18)\begin{equation} {\breve {\varPhi}} _{B} (\tilde{x}, y) = (1/T_{B}, {U_B}-c_r, 0, 0, T_{B})(\kern0.7pt y) + \sum_{m ={-} M_{W}}^{M_{W}} {{{\mathop{\phi}\limits^{{\smile}}}_{m0}}\left(\kern0.7pt y \right){\exp({\text{i} m \breve \alpha \tilde{x}}}}) +\cdots , \end{equation}

where $c_r= \omega _{0}/{\breve \alpha }$, $\tilde {x} = x - c_rt$ is the Galilean transformed coordinate and $M_W$ denotes the order of the Fourier-series truncation. The laminar base flow $\varPhi _{B}$ develops with a length scale much greater than the wavelength of the fundamental mode $2 {\rm \pi}/ {\breve \alpha }$, and so the non-parallelism of $\varPhi _B$ in the local region is neglected in the present analysis, rendering a periodic feature of $\breve \varPhi _B$ in the streamwise direction.

According to the Floquet theory, the periodic base flow supports the instability modes $\tilde {\phi }_W$ which can be expressed as

(2.19)\begin{equation} \left. \begin{gathered} {{\tilde{\phi}_W }(\tilde{x},y,z,t)} = \epsilon_W {\breve \phi _W}\left( {\tilde{x},y} \right){\exp{[\tilde{\sigma}( \tilde{x}+c_rt)+\text{i} \beta z+{\text{i}{\tilde{\sigma} _d} \breve \alpha \tilde{x}}]}} + \text{c.c.},\\ {\breve \phi _W}\left( {\tilde{x},y} \right) = \sum_{n ={-} {N_W}}^{{N_W}} {{{\hat \phi }_{W,n}}\left(\kern0.7pt y \right)\exp \left( { \text{i} n \breve \alpha \tilde{x}} \right)}, \end{gathered} \right\} \end{equation}

where $\tilde {\sigma }$ represents the growth rate, $\beta$ is the spanwise wavenumber, $\tilde {\sigma }_{d}$ is the detuning parameter, and $N_W$ is the order of the Fourier-series truncation and $\epsilon _W \ll 1$ measures the amplitude. For the fundamental resonance, we take $\tilde {\sigma } _d=0$. The component $\hat {\phi }_{W,0}$ denotes the streak component, and $\hat {\phi }_{W,n}$ with $n \ne 0$ represents the travelling mode. Substituting (2.18) and (2.19) into (2.7) with $O(\epsilon _W^{2})$ terms neglected, we arrive at a linear system,

(2.20)\begin{equation} ({{{{\boldsymbol M}_0}} + \tilde{\sigma} { {{\boldsymbol M}_1}} + {\tilde{\sigma}^2} { {{\boldsymbol M}_2}} }){ \breve \phi _W}\left( {\tilde{x},y} \right) = 0, \end{equation}

where

(2.21)\begin{equation} \left. \begin{aligned} { {{\boldsymbol M}_0}} & = \left( { {\boldsymbol{A}} + {\rm{i}} \beta { {{\boldsymbol{V}}_{ x z}}}} \right)\frac{\partial }{{\partial \tilde{x}}} + \left( {{\boldsymbol {B}} + {\rm{i}} \beta { {{\boldsymbol V}_{y z}}}} \right)\frac{\partial }{{\partial y}} + [ { {\boldsymbol {D}} + {\rm{i}} \beta {\boldsymbol { C}} + {{( {{\rm{i}} \beta } )}^2}{ { {{\boldsymbol V}_{ z z}}}}} ]\\ & \quad +{{{\boldsymbol V}_{xx}}}\frac{{{\partial ^2}}}{{\partial {{\tilde{x}}^2}}} + { {{\boldsymbol V}_{yy}}}\frac{{{\partial ^2}}}{{\partial {y^2}}} + { {{\boldsymbol V}_{xy}}}\frac{{{\partial ^2}}}{{\partial \tilde{x} \partial y}},\\ { {{\boldsymbol M}_1}} & = \left( { {\boldsymbol {A}} + \text{i} \beta { {{\boldsymbol V}_{ x z}}}} \right) + 2 { {{\boldsymbol V}_{xx}}}\frac{\partial }{{\partial \tilde{x}}} + {{{\boldsymbol V}_{xy}}}\frac{\partial }{{\partial y}} + c {\boldsymbol {G}} ,\\ { {{\boldsymbol M}_2}} & = { {{\boldsymbol V}_{xx}}}. \end{aligned} \right\} \end{equation}

The wall-normal boundary conditions read

(2.22)\begin{equation} \left. \begin{gathered} {\breve u}_{W}= {\breve v}_{W}={\breve w}_{W}={\breve T}_{W}=0, \quad \text{at} \ y=0, \\ ({\breve \rho}_{W},{\breve u}_{W}, {\breve v}_{W},{\breve w}_{W},{\breve T}_{W}) \to 0, \quad \text{as} \ y \to \infty. \end{gathered} \right\} \end{equation}

The linear system (2.20) with the homogeneous boundary conditions forms an eigenvalue problem with the growth rate $\tilde {\sigma }$ being the eigenvalue. In our paper, we choose $(M_W,N_W)=(3,6)$, which has been confirmed to be of sufficient accuracy. Such an analysis has also been used in the study of the SI of 2-D Mack second modes in hypersonic boundary layers (Chen et al. 2017; Xu et al. 2020), and our code validation and discretisation method are provided by Song, Zhao & Dong (2023).

3. Demonstration of the fundamental resonance regime by NPSE calculations

3.1. Case studies

For demonstration of the FR, we select a wind-tunnel condition of Maslov et al. (2001), for which the Mach number and temperature of the oncoming stream are 5.92 and 48.69 K, respectively. Such an oncoming condition was also used by Dong et al. (2020) and Dong & Zhao (2021). Two wall temperatures as listed in table 1 are selected, which are equal to and a half of the adiabatic wall temperature $T_{ad}$, where $T_{ad}$ is estimated by an empirical formula in White (2006, p. 512),

(3.1)\begin{equation} T_{ad}=1+\sqrt{Pr}(\gamma-1) M^{2}/2. \end{equation}

The nominal boundary-layer thickness for each case is also listed in the table.

Table 1. Parameters characterising the flow condition.

3.2. Base flow and its linear instability

The base-flow profiles of $U_B$ and $T_B$ at $x=0$ for the two cases are shown in figure 3. As the wall temperature decreases, the boundary-layer thickness is reduced and the shear rates of $U_B$ and $T_B$ at the wall increase. Solving the O-S equations numerically based on these base-flow profiles, we obtain the dependence of the growth rates $-\alpha _i$ of 2-D Mack modes on the frequency $\omega$ for the two cases, as shown in figures 4(a) and 4(b). Two distinguished unstable zones appear for each case, which are marked by the Mack first and second modes (Mack 1987). The second mode is more unstable than the first mode, and decrease of the wall temperature leads to an enhancement of the second mode and suppression of the first mode overall. In each panel, we show the results for three Reynolds numbers, namely, $R=2 \times 10^{3}$, $1\times 10^{4}$ and $\infty$ (for this case, the O-S equations reduce to the Rayleigh equations, which will be shown in § 4.1). Overall, increase of the Reynolds number leads to a greater growth rate, indicating the inviscid nature of the 2-D Mack modes.

Figure 3. (a) Streamwise velocity and (b) temperature of the compressible Blasius solution at $x=0$ forCases A and B. The red and black horizontal lines denote the nominal boundary-layer thicknesses for the two cases, respectively.

Figure 4. Dependence on the frequency $\omega$ of the growth rate $-\alpha _i$ of 2-D modes for (a) Case A and(b) Case B.

The accumulated amplitude of each Fourier mode can be quantified by an $N$ factor according to the LST, defined by

(3.2)\begin{equation} N(x) =\exp\left[ \int_{{0}}^x { - {\alpha _i} (\bar x) \,{\rm d}\kern0.06em \bar x}\right]. \end{equation}

In figure 5, we plot the streamwise evolution of the $N$ factors of 2-D second modes with representative frequencies in the second-mode frequency band for $R=2 \times 10^{3}$. It is observed that for Cases A and B, the frequencies of the most amplified second modes from $x=0$ to $1200$ are $\omega =0.11$ and $\omega =0.125$, respectively, and they are selected as the fundamental frequencies $\omega _0$ in the following NPSE calculations.

Figure 5. Streamwise evolution of the $N$ factors of 2-D second modes for (a) Case A and (b) Case B at $R=2 \times 10^{3}$.

3.3. Calculations of the fundamental resonance

For each case, we calculate the nonlinear evolution of the initial perturbations (2.8) using the NPSE approach, until the calculation blows up, indicating the emergence of the transition onset in a short distance downstream (Dong, Zhang & Zhou 2008). The parameters for Cases A and B are summarised in table 2. For Case A, the contours of the instantaneous velocity $u$ in the $x\unicode{x2013} z$ plane at two wall-normal positions are shown in figures 6(a) and 6(b). For $x<930$, the perturbation field is dominated by planar waves; however, 3-D structures appear at further downstream locations $(x>930)$ and grow with a high rate. Panels (c,d) show the time-averaged streamwise velocity at the same wall-normal positions for comparison, and the low- and high-speed streaks are observed evidently in the late nonlinear phase. In panel (e), we plot the contours of the time-averaged streamwise velocity in the $y$–$z$ plane at $x=1360$, where the spanwise localised blue structure indicates the low-speed streaks. The streamlines show the counter-rotating roll structures of the streamwise vorticities, which push the near-wall fluids upward, showing a lift-up mechanism for the formation of the low-speed streaks.

Table 2. Parameters for case studies in § 3.3.

Figure 6. Contours of the velocity field obtained by the NPSE for Case A: (a,b) instantaneous velocity $u$ in the $x$–$z$ planes at $y=5\,\% \delta _{99}$ and $y=40\,\% \delta _{99}$, respectively; (c,d) time-averaged velocity in the $x$–$z$ planes at $y=5\,\% \delta _{99}$ and $y=40\,\% \delta _{99}$, respectively; (e) time-averaged velocity in the $y$–$z$ plane at $x=1360$, where the arrayed curves show the streamlines.

The amplitude of each Fourier component in the physical space can be expressed as

(3.3)\begin{equation} {\tilde{u}_{\max }^{(m,n)}(x)} = \left\{ \begin{array}{@{}ll} {\max _y}\left| {{{\mathop{u}\limits^{{\smile}}}_{mn}}\left( {x,y} \right) + {\rm c.c.}} \right|, & {{ m \ne 0 \lor n \ne 0}},\\ {\max _y}\left| {{{\mathop{u}\limits^{{\smile}}}_{mn}}\left( {x,y} \right) } \right|, & {{ m=0 \land n=0}}. \end{array} \right. \end{equation}

In figure 7(a), we plot the streamwise evolution of ${\tilde {u}}_{max}^{(m,n)}$, shown by the solid lines. The amplitude of the fundamental mode $(1,0)$ agrees well with the linear result predicted by LPSE (shown by the red circles) until $x \approx 1100$, after which it saturates due to the nonlinearity. The oblique waves $(1,\pm 1)$ grow with exactly the same rate, and so only the curve for ${\tilde {u}}_{max}^{(1,1)}$ is plotted. It agrees with the linear prediction shown by the black circles until $x \approx 700$, after which a drastic amplification is observed. Although the other Fourier components are not introduced as initial perturbations, they are excited due to the mutual interaction of the introduced modes. The mean-flow distortion (MFD) $(0,0)$ and the harmonic mode $(2,0)$ are driven by the self-interaction of mode $(1,0)$, and therefore, their growth rates are almost twice that of $(1,0)$ in most of the computational domain, as confirmed by comparison with the blue crosses. For $x > 700$, the streak component $(0,1)$ and the high-order harmonics $(2,1)$ grow at almost the same rate as $(1,1)$, but their amplitudes differ by a remarkable amount; the amplitude of the streak mode $(0,1)$ is the greatest among the three. When the streak mode $(0,1)$ overwhelms the fundamental mode $(1,0)$ and becomes the dominant perturbation, the calculation blows up, indicating that the transition to turbulence is not far.

Figure 7. Streamwise evolution of the NPSE results for Case A: (a) amplitude of each Fourier component $\hat u_{max}^{(m,n)}$; (b) perturbation energy of each Fourier component $E_{mn}$; (c) energy growth rate of each Fourier component $\sigma _E$; (d) coefficient of surface friction, where the $C_f$ curve for turbulence is given by the empirical formula of White (2006, p. 553).

Alternatively, one can trace the evolution of the perturbation energy of each Fourier component, which is defined as (Chu 1965)

(3.4a,b)\begin{align} E_{mn}=\int_{0}^{\infty}{\mathop{\phi}\limits^{{\smile}}}_{mn}^{{\dagger}} M{{\mathop{\phi}^{{\smile}}}_{mn}} \,{{\rm d} y}, \quad M={\rm diag}\left(\frac{T_B}{ \gamma M^2 \rho_B}, \rho_B,\rho_B,\rho_B,\frac{\rho_B}{ \gamma (\gamma - 1) M^{2} T_B } \right) , \end{align}

where the superscript ${\dagger}$ denotes the complex conjugate with respect to its argument. The streamwise evolution of $E_{mn}$ for each Fourier component is shown in figure 7(b), and overall the same feature as in panel (a) is observed. This is quite predictable, because $\tilde {u}$ and $\tilde {T}$ are the dominant components in $\tilde {\phi }$ with similar evolution trend, and the evolution of the perturbation energy should agree overall with that of each dominant component. Figure 7(c) further plots the evolution of the growth rate of the perturbation energy, defined as

(3.5)\begin{equation} \sigma_E^{(m,n)} =\frac{1}{2E_{mn}} \frac{ {\rm d} E_{mn}}{{\rm d}\kern0.06em x}. \end{equation}

Remarkably, in the interval of $x \in [900,1200]$, as highlighted by the dashed box, the growth rates of $(0,1)$, $(1,1)$ and $(2,1)$ are almost identical, and much greater than that of the fundamental mode in the linear phase. This is a representative feature of the FR, as also observed by Sivasubramanian & Fasel (2015), Chen et al. (2017) and Hader & Fasel (2019). Such a high growth rate was explained by the SIA as introduced in § 2.3.3. In the SIA, the base flow is regarded as a superposition of the time- and spanwise-averaged mean flow and the quasi-saturated 2-D fundamental mode, together with its high-order harmonics, and the perturbation fields, including the streak mode, the 3-D travelling mode and higher-order harmonics with the same spanwise wavenumber, are governed by a linear eigenvalue system (2.20), with the growth rate $\tilde {\sigma }$ appearing as the eigenvalue.

In figure 7(d), the streamwise evolution of the coefficients of the skin friction,

(3.6)\begin{equation} C_f = \left( \frac{{2\underline{\mu}}}{R}\frac{{\partial \underline{u}}}{{\partial y}}\right)_{y=0}, \end{equation}

are plotted, where $\underline {u}$ and $\underline {\mu }$ represent the temporal and spanwise average of the streamwise velocity and the dynamic viscosity. The $C_f$ curve obtained by the NPSE calculation decreases with $x$ gradually at the beginning, agreeing with the unperturbed laminar-flow state, but it starts to deviate from the laminar state at $x \approx 500$, indicating a moderate MFD appearing there. The $C_f$ curve reaches its first peak at $x \approx 1000$, followed by a plateau until $x \approx 1300$, after which it shows another increase until the blowup position. The double-increase phenomenon is typical for the fundamental resonance regime, as also reported by previous works (Chen et al. 2017; Hader & Fasel 2019). The first increase is associated with the strong MFD induced by the finite-amplitude fundamental mode, and the following plateau agrees with the region of FR. Due to the FR, the streak mode becomes the dominant perturbation in the late phase, which, together with the travelling modes, may drive another type of secondary instability to support the growth of the high-frequency perturbations. Since these secondary instability modes amplify with high growth rates, which produce sufficient Reynolds stress to cause the rapid distortion of the mean flow, the parabolised assumption in the NPSE approach ceases to be valid, leading to the blowup of the NPSE calculation eventually.

In figure 8, we compare the growth rates of modes $(1,1)$ and $(0,1)$, $\sigma _E^{(1,1)}$ and $\sigma _E^{(0,1)}$, with that predicted by SIA, $\tilde {\sigma }$. In the interval $x \in [900,1200]$, the three curves agree with each other. The perturbation profiles of $(0,1)$ and $(1,1)$ obtained by the two approaches also agree well, as shown in figure 9. In fact, $\hat {u}_{01}$, $\hat {v}_{01}$ and $\hat {T}_{01}$ are real and $\hat {w}_{01}$ is pure imaginary, which will be discussed in § 4.2.1. Observations in figure 9 also indicate that the amplitude of the streamwise velocity of the streak mode $\hat {u}_{01}$ is much greater than that of the 3-D travelling mode $\hat {u}_{11}$, agreeing with the amplitude evolution in figure 7(a). For the profiles of the temperature perturbation, the streak mode in the bulk of the boundary layer $\hat {T}_{01}$ is also much greater than the 3-D travelling mode $\hat {T}_{11}$, but it is much weaker in the near-wall region. Although the SIA can predict quantitatively the growth rates and profiles of both the streak mode and 3-D travelling modes in the interval $[900,1200]$, the underlining mechanism determining the dominant role of the streak mode in the bulk region is not obvious. To answer these questions, a more in-depth analysis is required as will be introduced in the following section. Numerical results for Case B show the same feature, which will be illustrated in detail in § 5.

Figure 8. Comparison of the growth rate $\sigma _E$ obtained by the NPSE calculation with the SIA prediction $\tilde {\sigma }$.

Figure 9. Comparison of the perturbation profiles for $(0,1)$ and $(1,1)$ obtained by the SIA and NPSE approaches for Case A at $x=1000$. The profiles are normalised by the maximum of ${\hat v}_{01}$.

4. Asymptotic analysis for the principle of the fundamental resonance

4.1. Flow decomposition and the 2-D fundamental mode

To reveal the principle mechanism of the FR, we perform a weakly nonlinear analysis based on the high-Reynolds-number asymptotic technique. In the weakly nonlinear phase, the perturbation field $\tilde {\phi }=(\tilde {\rho },\tilde {u}, \tilde {v}, \tilde {w},\tilde {T})$ defined in (2.6) includes a set of harmonic perturbations,

(4.1)\begin{equation} \tilde{\phi}(x,y,z,t)=\bar{\epsilon}_{00}\tilde{\phi}_{00}+ \bar{\epsilon}_{10}\tilde{\phi}_{10} +\bar{\epsilon}_{11}\tilde{\phi}_{11}+\bar{\epsilon}_{1-1}\tilde{\phi}_{1-1}+\bar{\epsilon}_{01}\tilde{\phi}_{01}+\text{c.c.} +{\rm h.o.t.}, \end{equation}

where $\bar \epsilon _{mn} \ll 1$ denotes the amplitude of each component in the nonlinear phase, $\textrm {h.o.t.}$ denotes the high-order terms, and the subscripts $00$, 10, 11 ( and 1–1) and 01 denote the MFD, the fundamental mode, the 3-D travelling modes and the streak mode, respectively. According to the numerical results in § 3.3, we know that the amplitude of the fundamental mode $\bar {\epsilon }_{10}$ is much greater than those of the 3-D travelling modes and the streak mode, namely,

(4.2a,b)\begin{equation} \bar{\epsilon}_{10}\gg \bar{\epsilon}_{11},\quad \bar{\epsilon}_{10}\gg \bar{\epsilon}_{01}. \end{equation}

In the following analysis, we will focus on the region $x>900$, for which the fundamental mode $\tilde {\phi }_{10}$ evolves either in the linear phase or in the nonlinear phase, but the streak mode and the 3-D travelling modes undergo drastic amplification, as shown in figure 7(a). The base flow for the nonlinear analysis is chosen as the time- and spanwise-averaged mean flow, which includes the Blasius solution and the MFD,

(4.3)\begin{equation} (\bar U,\bar T)(x,y)=( U_B, T_B)(x,y)+\bar{\epsilon}_{00}(\tilde{u}_{00},\tilde{T}_{00})(x,y). \end{equation}

In the early nonlinear phase, the MFD is mainly driven by the fundamental mode, which also acts back on the fundamental mode to lead to its saturation eventually. Since the streamwise length scale of the mean flow is much greater than the Mack-mode wavelength, the non-parallelism of the base flow is negligible in the following analysis.

From the linear stability analysis based on the parallel mean flow (4.3) at a chosen streamwise location $x$, we find that the growth rate of the fundamental Mack mode is much smaller than its wavenumber. We can express the perturbation profiles $\tilde {\phi }_{10}$ in terms of

(4.4)\begin{equation} \tilde{\phi}_{10}(x,y,t)=\hat \phi_{10}(\kern0.7pt y;x) \exp\left({{\rm i}\left(\int^x_{0}\alpha \,{{\rm d}\kern0.06em x}-\omega_0 t\right)}\right), \end{equation}

where the streamwise wavenumber $\alpha$ is almost real and $O(1)$, the frequency $\omega _0$ is also taken to be $O(1)$ and $\hat \phi _{10}$ is the eigenfunction of the fundamental mode. Asymptotic analyses, as done by Dong et al. (2020) and Dong & Zhao (2021), showed that the Mack mode shows a double-deck structure in the high-$R$ approximation, namely, a main layer where $y=O(1)$ and a viscous Stokes layer where $y=O(R^{-1/2})$.

The eigenfunction of the fundamental mode $\hat {\varphi }_{10} = (\hat {v}_{10},\hat {p}_{10})$ satisfies the Rayleigh equation in the main layer based on the mean flow $(\bar U, \bar T)$,

(4.5)\begin{equation} {\mathcal{L}}_{R} \hat{\phi}_{10} \equiv (d_y - {\mathcal{H}}_{0} ) \hat{\varphi}_{10}=0 , \end{equation}

where

(4.6)\begin{equation} {\mathcal{H}}_{0}= \left( \begin{array}{@{}cc@{}} {{\bar U_{y}}}/({{\bar U - c}}) & {{ - {M^2}S_0 - {\alpha ^2}\bar T}}/{{{S_0}}}\\ - {{{S_0}}}/{{\bar T}} & 0 \end{array} \right) \end{equation}

with $S_0 = \textrm {i} \alpha (\bar U -c)$ and $c \equiv \omega /\alpha$. The boundary conditions read

(4.7a,b)\begin{equation} \hat v_{10}(0)=0;\quad \hat p_{10}(\infty)\rightarrow 0. \end{equation}

Such a linear, homogeneous system leads to an eigenvalue problem. For a spatial mode, the frequency $\omega$ is given to be real, and following the numerical method as used by Dong et al. (2020), Dong & Zhao (2021) and Zhao & Dong (2022), we can calculate the eigenvalue $\alpha =\alpha _r+\textrm {i} \alpha _i$ $(|\alpha _i| \ll |\alpha _r|)$ and the eigenfunctions. If, for a particular $\omega$, the Mack mode is neutral, i.e. $\alpha _i=0$, then the system (4.5) becomes singular, and a critical layer around the location where $\bar {U}=\omega /\alpha$ appears. For a linear critical layer, the leading-order balance in this layer is between the inertial and viscous terms, and in the numerical process, the solution can be obtained by detouring the integrating path around the critical point, as illustrated by Schmid & Henningson (2001, pp. 44–45); an example can be found from Dong et al. (2020). Other perturbation quantities such as $\hat {u}_{10}$ and $\hat {T}_{10}$ can be obtained by

(4.8a,b)\begin{equation} \hat{u}_{10}={-}\frac{\bar U_{y} \hat{v}_{10}+ {\rm i} \alpha \hat{p}_{10}}{S_{0}},\quad \hat{T}_{10}={-}\frac{\bar T_{y} \hat{v}_{10}}{S_{0}}+(\gamma-1)M^{2}\bar T \hat{p}_{10}. \end{equation}

Additionally, we take $\textrm {max}_{y} |\hat {u}_{10}(\kern0.7pt y)|=1$ for normalisation.

Because the Rayleigh solution does not satisfy the no-slip condition at the wall, a viscous Stokes layer needs to be taken into account, whose solutions can be found from Dong et al. (2020) and so are not presented here.

4.2. The 3-D travelling modes and the streak mode

It is seen from figure 7(a) that after the fundamental 2-D mode reaches a finite amplitude, both the 3-D travelling modes $(1,\pm 1)$ and the streak mode $(0,1)$ amplify with the same rate, which is much greater than that of the fundamental mode, showing an FR phenomenon. To reveal this mechanism, we probe the evolution of the perturbations $\tilde {\phi }_{1\pm 1}$ and $\tilde {\phi }_{01}$, which are expressed as

(4.9a,b)\begin{gather} \begin{array}{c} \tilde{\phi}_{11}=\hat\phi_{11}(\kern0.7pt y) \exp({{\rm i}(\alpha_r x+\beta_0 z-\omega_0 t)+\sigma x}),\\ \tilde{\phi}_{1-1}=\hat\phi_{1-1}(\kern0.7pt y) \exp({{\rm i}(\alpha_r x-\beta_0 z-\omega_0 t)+\sigma x}),\end{array} \end{gather}

(4.10)\begin{gather} \tilde{\phi}_{01}=\hat\phi_{01}(\kern0.7pt y)\exp({{\rm i}\beta_0 z+\sigma x}), \end{gather}

where $\sigma$ denotes their common growth rate. Due to the fundamental resonance, the streamwise wavenumber $\alpha _r$ of modes $(1,\pm 1)$ is the same as that of the fundamental 2-D mode, and all these modes grow with the same rate $\sigma$. Although the growth rate $\sigma$ could be much greater than that of the fundamental mode, its magnitude is still much smaller than unity, i.e.

(4.11)\begin{equation} |\alpha_i|\ll \sigma\ll 1. \end{equation}

4.2.1. Leading-order asymptotic solutions in the main layer

In the main layer where $y\sim 1$, balancing the governing equations and taking into account the scalings of $\alpha$, $\beta$ and $\omega$, we obtain that all the quantities of the 3-D travelling mode, $\hat {\boldsymbol {u}}_{1 \pm 1}$, $\hat \rho _{1 \pm 1}$, $\hat T_{1 \pm 1}$, $\hat p_{1 \pm 1}$ are of the same order, but those of the streak mode are not. The streamwise wavenumber of the streak mode is zero, and its streamwise growth rate $\sigma$ is small. Thus, the streamwise derivative of $\hat \phi _{01}$ with respect to $x$ is only $O(\sigma \hat \phi _{01})$, much smaller than its derivatives with respect to $y$ and $z$. Then, from the balance of the continuity equation, we obtain that the streamwise perturbation $\hat u_{01}$ must be much greater than $\hat v_{01}$ and $\hat w_{01}$. For convenience, we let $\hat u_{01}\sim 1$ (because its magnitude is measured by $\bar \epsilon _{01}$), then balance of the continuity equation leads to

(4.12)\begin{equation} \hat v_{01}\sim\hat w_{01}\sim \sigma. \end{equation}

Here, $\tilde {\phi }_{11}$ is driven by the nonlinear interaction of $\tilde {\phi }_{01}$ and $\tilde {\phi }_{10}$, therefore, we have

(4.13)\begin{equation} \bar{\epsilon}_{11}=\bar{\epsilon}_{10} \bar{\epsilon}_{01}. \end{equation}

In the spanwise momentum equation of the streak mode, the inertial term $uw_x\sim \bar {\epsilon }_{01}\sigma \bar U \hat w_{01}$, the pressure-gradient term $p_z\sim \bar {\epsilon }_{01}\textrm {i}\beta \hat p_{01}$ and the nonlinear terms $O(\bar {\epsilon }_{10} \bar {\epsilon }_{11})$ should balance, namely,

(4.14)\begin{equation} \bar\epsilon_{01}\sigma^2\sim\bar\epsilon_{01}\hat p_{01}\sim \bar\epsilon_{10}\bar\epsilon_{11}, \end{equation}

leading to

(4.15a,b)\begin{equation} \sigma\sim \bar{\epsilon}_{10},\quad \hat p_{01}\sim \sigma^{2}\sim \bar{\epsilon}_{10}^2. \end{equation}

In the balance of the streamwise momentum equation, the nonlinear terms $O(\bar \epsilon _{10}\bar \epsilon _{11})$ are much smaller than the inertia term $O(\bar \epsilon _{01}\sigma )$, and so they do not appear in the leading-order balance. For convenience, we introduce

(4.16)\begin{equation} \bar \sigma=\bar{\epsilon}_{10}^{ {-}1}{\sigma}=O(1). \end{equation}

It should be noted that one may be dissatisfied with (4.15a) if the NPSE results are re-examined, because figure 7(a) shows that in the FR region, $\bar \epsilon _{10} \sim 0.1$, but figure 8 indicates that $\sigma \sim 0.01$. Actually, the scaling relation (4.15a) is from a priori analysis from the properties of the governing equations, instead of an observation from the numerical results at a finite Reynolds number. The purpose to prescribe this scaling relation is to construct a consistent mathematical description to explain the physics of FR, which should agree with the reality as long as the Reynolds number is sufficiently high. Therefore, the ratio of $\sigma$ and $\bar \epsilon _{01}$ is expected to approach $O(1)$ as $R$ increases, which will be displayed in figure 22.

Likewise, balancing the energy equation and the equation of the state for $\hat \phi _{01}$, we obtain that $\hat T_{01}\sim \hat \rho _{01}\sim 1$. Thus, we introduce

(4.17)\begin{equation} (\hat u_{01},\hat v_{01},\hat w_{01},\hat p_{01},\hat \rho_{01},\hat T_{01})=(\breve u_{01},{\bar{\epsilon}_{10}}\breve v_{01},\bar{\epsilon}_{10}\breve w_{01},\bar{\epsilon}_{10}^{2}\breve p_{01},\breve \rho_{01},\breve T_{01})+\cdots. \end{equation}

The physical quantities for the 3-D travelling waves $\hat \phi _{1\pm 1}$ are all $O(1)$. For convenience, both $\hat {\phi }_{1+1}$ and $\hat {\phi }_{11}$ are used to denote the 3-D travelling mode $(1,1)$ in this paper.

Collecting the leading-order terms from the governing equations of the 3-D travelling modes $\tilde {\phi }_{1 \pm 1}$, we obtain

(4.18a) \begin{gather} \hat{S}_0 M^2\hat p_{1 \pm 1}+{\rm i}\alpha_r \hat u_{1 \pm 1}+\hat v_{1 \pm 1}' \pm {\rm i}\beta \hat w_{1 \pm 1}=F_{1 \pm 1,1} \nonumber\\ \equiv{-} {M^2}{\rm i}\alpha_r {\hat p_{10}} {\breve u}_{01} + \left( {\hat{S}_0{{\hat \rho }_{10}} - \frac{{\bar T_{y}}}{{{{\bar T}^2}}}{{\hat v}_{10}}} \right) {\breve T}_{01} - \bar T( {{\rm i}\alpha_r {{\hat u}_{10}} + {{\hat v'}_{10}}}) {\breve \rho}_{01} , \end{gather}

(4.18b) \begin{gather} \hat{S}_0 \hat u_{1 \pm 1}+\bar U_{y}\hat v_{1 \pm 1}+{\rm i}\alpha_r \bar T\hat p_{1 \pm 1}=F_{1 \pm 1,2} \equiv{-}\hat v_{10}\breve u_{01}'-{\rm i}\alpha_r \hat u_{10}\breve u_{01}+{\rm i}\alpha_r\bar T^2\hat p_{10}\breve \rho_{01}, \end{gather}

(4.18c) \begin{gather} \hat{S}_0 \hat v_{1 \pm 1}+\bar T\hat p_{1 \pm 1}'=F_{1 \pm 1,3}\equiv{-}{\rm i}\alpha_r \hat v_{10}\hat u_{01}-\hat{S}_0 \bar T\hat v_{10}\hat \rho_{01}, \end{gather}

(4.18d) \begin{gather} \hat{S}_0 \hat w_{1 \pm 1} \pm {\rm i}\beta\bar T\hat p_{1 \pm 1}=F_{1 \pm 1,4}\equiv 0, \end{gather}

(4.18e) \begin{gather} \hat{S}_0 {{\hat T }_{1 \pm 1}} + \bar T_{y}{{\hat v}_{1 \pm 1}} - \left( {\gamma - 1} \right){M^2}\bar T \hat{S}_0 {{\hat p}_{1 \pm 1}} =F_{1 \pm 1,5} \nonumber\\ \equiv{-} {{\hat v}_{10}} {\breve T '}_{01} - {\rm i}\alpha_r [ {{{\hat T }_{10}} - ( {\gamma - 1}){M^2}\bar T{{\hat p}_{10}}} ] {\breve u}_{01} - \bar T({\bar T_{y}{{\hat v}_{10}}-\hat{S}_0 {{\hat T }_{10}}}) \breve \rho_{01}, \end{gather}

(4.18f) \begin{gather} - \frac{{{{\hat T }_{1 \pm 1}}}}{{\bar T}} + \gamma {M^2}{\hat p_{1 \pm 1}} - \bar T{\hat \rho _{1 \pm 1}} = F_{1 \pm 1,6}\equiv {\hat \rho _{10}} {\breve T}_{01} + {\hat T _{10}} {\breve \rho}_{01}, \end{gather}

where $\hat {S}_0 = -\textrm {i}\omega +\textrm {i}\alpha _r \bar U$ and in what follows, the prime denotes the derivative with respect to its argument. Equation (4.18a) is obtained by eliminating $\hat \rho _{1\pm 1}$ and $\hat T_{1\pm 1}$ from the continuity equation, the energy equation, and the equation of the state. It is seen that a critical layer appears at a position where $\bar U=\omega /\alpha _r$, and we use the standard treatment for the linear critical layer as done by Schmid & Henningson (2001, pp. 44–45) to avoid the singularity.

If the profiles of the streak mode $\breve \phi _{01}$ are known, then (4.18) forms an inhomogeneous linear system. However, they are coupled with the unknown vector of the streak mode $\breve \phi _{01}$, rendering a triad resonance system.

The leading-order governing equations for the streak mode are

(4.19a)\begin{gather} \bar\sigma\breve u_{01}+\breve v_{01}'+{\rm i}\beta\breve w_{01}=0, \end{gather}

(4.19b)\begin{gather} \bar\sigma\bar U\breve u_{01}+\bar U_{y}\breve v_{01}=0, \end{gather}

(4.19c)\begin{gather} \bar\sigma\bar U\breve v_{01}+\bar T\breve p_{01}'= F_{01,3}\equiv{-}\hat v_{10}^{\dagger} (\hat v_{11}'-{\rm i}\alpha_r\hat u_{11}-\bar T \hat{S}_{0} \hat{ \rho}_{11})-\left(2(\hat v_{10}^{\dagger})'-\frac{\bar T_{y}}{\bar T}\hat v_{10}^{\dagger}\right)\hat v_{11}+{\rm c.c.}, \end{gather}

(4.19d)\begin{gather} \bar\sigma\bar U\breve w_{01}+{\rm i}\beta\bar T\breve p_{01}= F_{01,4}\equiv{-}(\hat v_{10}^{\dagger}\hat w_{11})'+\frac{\bar T_{y}}{\bar T}\hat v_{10}^{\dagger}\hat w_{11}-{\rm c.c.}, \end{gather}

(4.19e,f)\begin{gather} \bar\sigma\bar U\breve T_{01}+\bar T_{y}\breve v_{01}=0, \quad \breve T_{01}+\bar T^2\breve \rho_{01}=0. \end{gather}

From the symmetric feature of these equations, we know that $\breve u_{01}$, $\breve v_{01}$, $\breve T_{01}$ and $\breve p_{01}$ are real, and $\breve w_{01}$ is pure imaginary. However, the detouring method near the critical layer is used and thus breaks the symmetric feature slightly. However, this effect is very weak and, therefore, to the leading-order approximation, we also consider $\breve u_{01}$, $\breve v_{01}$, $\breve T_{01}$ and $\breve p_{01}$ to be real, and $\breve w_{01}$ to be pure imaginary.

Remarkably, it is seen from (4.19) that the transverse and lateral velocities of the streak mode are forced by the nonlinear interaction of the 2-D fundamental and 3-D travelling modes, while its streamwise velocity with a greater amplitude is driven by the linear lift-up mechanism. Such a mechanism is reminiscent of the stronger amplification of the streak mode in the oblique breakdown regime observed by Song et al. (2022).

The attenuation condition is imposed in the far field, and the analysis of the wall layer, as will be shown in § 4.3, implies that the non-penetration condition is imposed at the lower boundary,

(4.20a–c)\begin{gather} (\,\hat p_{11},\hat p_{1-1},\breve p_{01})\rightarrow 0,\quad\text{as}\ y\rightarrow \infty, \end{gather}

(4.20d–f)\begin{gather} \hat v_{11}(0)=\hat v_{1-1}(0)=\breve v_{01}(0)=0. \end{gather}

Combining (4.18) and (4.19), we obtain a six-order linear differential system,

(4.21)\begin{equation} \frac{{\rm d} \hat\varphi}{{{\rm d} y} }=\boldsymbol{A}\hat\varphi, \end{equation}

where $\hat \varphi =(\hat v_{11},\hat p_{11},\hat v_{1-1},\hat p_{1-1},\breve v_{01},\breve p_{01})$ and $\boldsymbol {A}$ can be deduced readily from (4.18) and (4.19). Such a homogeneous system with the homogeneous boundary conditions forms an eigenvalue system, with $\bar \sigma$ being the eigenvalue. The numerical approaches as used by Malik (1990) can be employed to solve this system.

Now we are interested in the near-wall behaviours of the perturbation fields for components $(0,1)$ and $(1,1)$. Although $\hat v_{11}$ and $\breve v_{01}$ satisfy the non-penetration condition at the wall, the perturbations of the streamwise velocity, spanwise velocity and temperature may be finite or even blow up, disagreeing with the no-slip and isothermal boundary conditions. From a scaling estimate, we find that as $y\to 0$,

(4.22)\begin{gather} \breve v_{01}\sim y\ln y,\quad (\breve u_{01},\breve T_{01})\sim \ln y,\quad \breve w_{01}\sim 1,\quad \breve p_{01}\sim 1, \end{gather}

(4.23)\begin{gather} \hat v_{11}\sim y\ln y,\quad(\hat u_{11},\hat T_{11})\sim \ln y, \quad\hat w_{11}\sim 1,\quad \hat p_{11}\sim 1. \end{gather}

In fact, the inhomogeneous forcing terms in (4.19d) in the vicinity of the wall are expanded as

(4.24)\begin{equation} F_{01,4}=\hat A_2+\hat A_3y+\cdots, \end{equation}

where

(4.25)\begin{gather} \hat A_2={-}2{\rm i}\mathrm{Im}[\lambda_v^{\dagger}\hat w_{11}(0)],\quad \hat A_3={-}2{\rm i}\mathrm{Im}\left[2\lambda_v^{\dagger}\hat w_{11}'(0)+\left(\bar\lambda_v^{\dagger}-\frac{\lambda_T\lambda_v^{\dagger}}{T_w}\right)\hat w_{11}(0)\right], \end{gather}

(4.26)\begin{gather} \hat w_{11}(0)=\frac{\beta T_w}{\omega}\hat p_{11}(0),\quad \hat w_{11}'(0)=\left[\frac{\alpha \beta\lambda T_w}{\omega^2}+\frac{\beta \lambda_T}{\omega}\right]\hat p_{11}(0), \end{gather}

with $\lambda _v=\hat v_{10,y}(0)$, $\bar \lambda _v=\hat v_{10,yy}(0)$, $\lambda =\bar U_{y}(0)$ and $\lambda _{T}=\bar T_{y}(0)$. Thus, applying $y\rightarrow 0$ to (4.19), we have

(4.27a)\begin{gather} \breve v_{01}'=\frac{\breve v_{01}}{y}+\frac{{\rm i}\beta}{\bar\sigma \lambda y}[{\rm i}\beta (T_w+\lambda_Ty)\breve p_{01}-\hat A_2]-\frac{{\rm i}\beta\hat A_3}{\bar\sigma \lambda}+\cdots, \end{gather}

(4.27b)\begin{gather} \breve p'_{01}={-}\frac{\bar\sigma\lambda y}{T_w}\breve v_{01}+\cdots. \end{gather}

The solutions of the transverse velocity and pressure of the streak mode are

(4.28)\begin{equation} \breve v_{01}\to\frac{{\rm i}\beta }{\bar\sigma\lambda}\left(-\hat A_3+\frac{\lambda_T\hat A_2}{T_w}\right)y\ln y+\cdots, \quad \breve p_{01}\rightarrow \frac{\hat A_2}{{\rm i}\beta T_w}+\cdots, \end{equation}

and the other perturbations behave like

(4.29)\begin{align} (\breve u_{01},\breve T_{01})\to\left(1,\frac{\lambda_T}{\lambda}\right)\frac{{\rm i}\beta }{\bar\sigma^2\lambda}\left(\hat A_3-\frac{\lambda_T\hat A_2}{T_w}\right)\ln y+\cdots,\quad \breve w_{01}\to \frac{1}{\lambda\bar\sigma}\left(\hat A_3-\frac{\lambda_T\hat A_2}{T_w}\right) +\cdots. \end{align}

Obviously, $\breve u_{01}$ and $\breve T_{01}$ are unbounded at the wall, which requires consideration of the viscosity in a thin layer, referred to as the viscous wall layer.

4.3. Viscous wall layer

Since the streak mode is singular at the wall, a viscous wall layer has to be taken into account. Balancing the inertial and viscous terms, we obtain the thickness of the wall layer, $y\sim (\bar {\epsilon }_{10}R)^{-1/3}$. For convenience, we introduce a local coordinate

(4.30)\begin{equation} Y=\epsilon^{{-}1} y=O(1),\quad \epsilon\equiv (\bar{\epsilon}_{10}R)^{{-}1/3}. \end{equation}

From the main-layer estimate (4.29), we obtain the magnitude of the streak mode in the wall layer,

(4.31)\begin{equation} (\breve u_{01},\breve v_{01},\hat p_{01},\breve T_{01},\breve \rho_{01})\sim (\ln \epsilon,\epsilon \ln \epsilon, 1,\ln\epsilon,\ln\epsilon). \end{equation}

However, to satisfy the wall-layer governing equations, the perturbation spanwise velocity must come to the leading-order continuity equation. Thus, we let

(4.32)\begin{equation} \breve w_{01}\sim\ln \epsilon, \end{equation}

which decays algebraically as $Y\to \infty$. In the spanwise momentum equation, the leading-order balance is between the pressure gradient and the inhomogeneous forcing, and the inertial term $\sigma \bar U\hat w_{01}$ appears only in the second order, which must balance with the second-order pressure gradient $\textrm {i} \beta \hat p_{01}^{(2)}$, where $\hat p_{01}^{(2)}$ denotes the second-order pressure perturbation. This balance leads to $\hat p_{01}^{(2)}\sim \bar \epsilon _{10}^2\epsilon \ln \epsilon$, and as $Y\to \infty$, the spanwise velocity perturbation $\hat w_{01}\to \bar \epsilon _{10}\ln \epsilon Y^{-1}$ or $\breve w_{01}\to \ln \epsilon Y^{-1}$.

Therefore, the streak-mode perturbation field is expanded as

(4.33a) \begin{gather} (\breve u_{01},\breve v_{01},\breve w_{01},\breve T_{01})=\ln\epsilon(\breve U_0,\epsilon\breve V_0,\breve W_0,\breve T_0)+(\breve U_1,\epsilon\breve V_1,\breve W_1,\breve T_1)+\cdots, \end{gather}

(4.33b)\begin{gather} \breve p_{01}=\breve P_0+\epsilon\ln\epsilon \breve P_1+\cdots, \end{gather}

where $\breve P_0=\hat A_2/(\textrm {i}\beta T_w)$. The presence of $\epsilon \ln \epsilon \breve P_1$ and $\epsilon \breve V_1$ induces an $O(\epsilon \ln \epsilon )$ correction to the main-layer solution.

The perturbation quantities of the corresponding 3-D travelling mode in the wall layer are expanded as

(4.34)\begin{equation} (\breve u_{11},\breve v_{11},\breve w_{11},\breve p_{11},\breve T_{11},\breve \rho_{11})=(\ln \epsilon\hat U_{0},\epsilon\ln \epsilon\hat V_{0},\hat W_{0},\hat P_{0},\ln \epsilon\hat T_{0},\ln \epsilon\hat R_{0})+\cdots. \end{equation}

The leading-order governing equations for the streak mode read

(4.35a)\begin{gather} \bar\sigma\breve U_0+\breve V_0'+{\rm i}\beta \breve W_0=0, \end{gather}

(4.35b,c)\begin{gather} \lambda\bar\sigma Y\breve U_0+\lambda \breve V_0-C_w\breve U_0''=0,\quad \breve P_{1}'=0, \end{gather}

(4.35d,e)\begin{gather} {\rm i}\beta T_w\breve P_0=\hat A_2,\quad \lambda\bar\sigma Y\breve W_0-C_w\breve W_0''+{\rm i}\beta T_w \breve P_1=0, \end{gather}

(4.35f)\begin{gather} \lambda \bar \sigma Y \breve T_0+\lambda_T\breve V_0-\frac{C_w}{Pr}\breve T_0''=0, \end{gather}

where $C_w=\mu _wT_w$. This is a linear homogeneous system, in which the nonlinear interaction of the fundamental mode and the 3-D travelling mode does not appear in the leading-order balance. Because the viscosity appears in the leading order of the streak-mode equations, the no-slip conditions are satisfied at the lower boundary,

(4.36)\begin{equation} \breve U_0(0)=\breve V_0(0)=\breve W_0(0)=\breve T_0(0)=0. \end{equation}

In the upper limit, the perturbation field must match the main-layer solutions, namely,

(4.37)\begin{equation} (\breve U_0,\breve V_0,\breve W_0,\breve T_0)\rightarrow \left(1,-\bar\sigma Y,O(1/Y),\frac{\lambda_T}{\lambda} \right)\frac{{\rm i}\beta}{\lambda\bar\sigma^2}\left(\hat A_3-\frac{\lambda_T\hat A_2}{T_w}\right)\quad\mbox{as }Y\rightarrow \infty. \end{equation}

From (4.35d,e), the matching condition and the no-slip condition ($\breve W_0(0)=0$), we find that

(4.38)\begin{equation} \breve P_0=\frac{\hat A_2}{{\rm i}\beta T_w},\quad \breve W_0={-}\frac{{\rm i} {\rm \pi}\beta T_w\breve P_1}{C_w^{1/3}(\lambda\bar\sigma)^{2/3}}\left(\mbox{Gi}(\eta)-\frac{\mbox{Gi}(0)}{\mbox{Ai}(0)}\mbox{Ai}(\eta)\right), \end{equation}

where $\eta =(\lambda \bar \sigma /C_w)^{1/3}Y$, $\mbox {Ai}$ and $\mbox {Gi}$ are the Airy's functions of the first kind and the Scorer's function, respectively; see Abramowitz & Stegun (1964, p. 448).

Equating (4.35a,b,e) and eliminating $\breve U_0$, $\breve W_0$ and $\breve P_0$, we obtain

(4.39)\begin{equation} C_w\breve V_0^{(4)}-\lambda\bar\sigma Y\breve V_0''=0, \end{equation}

whose solution reads

(4.40)\begin{equation} \breve V_0=\bar C\int_0^{\eta}\int_0^{\bar\eta}\mbox{Ai}(\hat\eta)\,{\rm d}\hat\eta \,{\rm d}\bar\eta, \end{equation}

where $\bar C$ is a constant to be determined later. The upper limit of $\breve V_0$ is

(4.41)\begin{equation} \breve V_0\to \left[\frac{1}{3}(\lambda\bar\sigma/C_w)^{1/3}Y-0.2588\right]\bar C,\quad\mbox{as }Y\rightarrow \infty. \end{equation}

Comparing with the matching condition (4.37), we obtain

(4.42)\begin{equation} {\bar C}={-}\frac{3{\rm i}\beta C_w^{1/3}}{(\lambda\bar\sigma)^{4/3}}\left(\hat A_3-\frac{\lambda_T\hat A_2}{T_w}\right). \end{equation}

The implication is that in the main layer, the second-order perturbation is driven by an outflux $\epsilon \ln \epsilon V_\infty$ with

(4.43)\begin{equation} V_\infty={-}0.2588\bar C=\frac{0.7764{\rm i}\beta C_w^{1/3}}{(\lambda\bar\sigma)^{4/3}}\left(\hat A_3-\frac{\lambda_T\hat A_2}{T_w}\right). \end{equation}

Thus, to predict the growth rate $\bar \sigma$ more accurately, an improved boundary condition will be introduced in § 4.4.

Substituting (4.40) into the continuity equation (4.35a), we obtain

(4.44)\begin{equation} \breve U_0={-}\frac{1}{\bar\sigma}\left((\lambda\bar\sigma/C_w)^{1/3}\bar C\int_0^\eta\mbox{Ai}(\bar\eta)d\bar\eta+{\rm i}\beta\breve W_0\right). \end{equation}

The energy equation (4.35e) leads to

(4.45)\begin{equation} \breve T_0=\frac{Pr^{1/3}\lambda_T{\rm \pi}}{C_w^{1/3}\lambda\bar\sigma} \left[-\mbox{Ai}\int_0^{\xi}\breve V_0\mbox{Bi}( \bar \xi)d \bar \xi+\mbox{Bi}\int_0^\xi \breve V_0\mbox{Ai}(\xi)\,{\rm d}\xi\right], \end{equation}

where $\xi =(\lambda \bar \sigma Pr/C_w)^{1/3} Y$. Applying (4.35b) at $Y=0$, we obtain $\breve U_0''(0)=0$, which leads to

(4.46)\begin{equation} \breve P_1={-}\frac{3{\rm i} \mbox{Ai}'(0)C_w^{1/3}}{T_w\beta(\lambda\bar\sigma)^{1/3}}\left(\hat A_3-\frac{\lambda_T\hat A_2}{T_w}\right). \end{equation}

The governing equations for the 3-D travelling mode are

(4.47a)$$\begin{gather} {\rm i}\alpha\hat U_{0}+\hat V_{0}'={-}{\rm i}\alpha M^2\hat p_{w0}\breve U_{0}-{\rm i}\omega \hat\rho_{w0}\breve T_0+({\rm i}\alpha \hat u_{w0}+\lambda_v)T_w^{{-}1}\breve T_0 , \end{gather}$$

(4.47b)$$\begin{gather}-{\rm i}\omega\hat U_{0}={-}{\rm i}\alpha\hat u_{w0}\breve U_{0}-\lambda_v Y\breve U_0'-{\rm i}\omega T_w^{{-}1}\hat u_{w0}\breve T_0, \end{gather}$$

(4.47c)$$\begin{gather}{\rm i}\omega \hat W_0+{\rm i}\beta T_w\hat P_0=0, \end{gather}$$

(4.47d)$$\begin{gather}-{\rm i}\omega \hat T_{0}={-}{\rm i}\alpha [\hat \theta_{w0}-(\gamma-1)M^2T_w\hat p_{w0}]\breve U_0-\lambda_vY\breve T_0'-{\rm i}\omega T_w^{{-}1}\hat \theta_{w0}\breve T_0, \end{gather}$$

where $(\,\hat p_{w0},\hat u_{w0},\hat \theta _{w0},\hat \rho _{w0})=(\,\hat {p}_{10}(0),\hat u_{10}(0),\hat T_{10}(0),\hat \rho _{10}(0))$. It is seen that $\breve P_0$, $\breve V_0$, $\breve W_0$ and $\breve P_1$ do not appear in the leading-order balance. From the Rayleigh equation (4.5), we know that

(4.48)\begin{equation} \left(\hat u_{w0},\lambda_v,\hat \rho_{w0},\hat \theta_{w0}\right)=\left(\frac{\alpha T_w}{\omega},\left({\rm i}\omega M^2+\frac{\alpha^2T_w}{{\rm i}\omega}\right),\frac{M^2}{T_w},(\gamma-1)M^2T_w\right) \hat{p}_{w0}. \end{equation}

The viscous effect of the 3-D travelling mode is secondary in the wall layer, and so only the non-penetration condition is satisfied at the lower boundary,

(4.49)\begin{equation} \hat V_0(0)=0. \end{equation}

Solving the above system, we obtain the solutions for the streamwise velocity and temperature of the 3-D travelling mode,

(4.50)\begin{gather} \hat U_{0}=\left[ \frac{\alpha^2 T_w}{\omega^2}\breve U_{0}+\left( M^2-\frac{\alpha^2T_w}{\omega^2}\right) Y\breve U_0'+\frac{\alpha}{\omega}\breve T_0 \right] \hat{p}_{w0}, \end{gather}

(4.51)\begin{gather} \hat T_{0}=\left[ \left(M^2-\frac{\alpha^2T_w}{\omega^2}\right)Y\breve T_0'+ (\gamma-1)M^2\breve T_0 \right] \hat{p}_{w0}. \end{gather}

It is seen that the no-slip and isothermal conditions are satisfied automatically, i.e. $\hat U_0(0)=\hat T_0(0)=0$. Integrating (4.47a), we obtain the solution for the transverse velocity,

(4.52)\begin{equation} \hat V_0=\int_0^Y -{\rm i} \alpha \hat{p}_{w0} \left[ \left(M^2 +\frac{\alpha^2}{\omega^2} \right)\breve U_0 +\left(M^2 -\frac{\alpha^2}{\omega^2} \right) \bar Y \breve U_0' +\frac{\alpha}{\omega} \breve T_0 \right] \,{\rm d} \bar Y. \end{equation}

In the upper limit, we can estimate, by dropping the $O(Y)$ unbounded part, the outflux to the main layer,

(4.53)\begin{align} V_{1\infty}\equiv \frac{\alpha \beta \hat{p}_{w0} C_{w}^{1/3} }{{\bar \sigma}^{7/3} \lambda^{4/3} } \left[ C_1 \left(M^2+\frac{\alpha^2}{\omega^2} \right)+ C_2 \left(M^2-\frac{\alpha^2}{\omega^2} \right) +C_3 \frac{\lambda_T \alpha }{\lambda \omega} \right]\left(\hat A_3-\frac{\lambda_T\hat A_2}{T_w}\right). \end{align}

We obtain from numerical calculations that $C_1=-0.8869$, $C_2=0.1139$ and $C_3=-1.0221$.

Interestingly, from (4.51), we find that the second term of $\hat T_{0}$ appears as $(\gamma -1)M^{2}\breve T_{0}$. Although we take $M=O(1)$ in this paper, the factor $M^2$ may be numerically large for a hypersonic case. Thus, for a finite-$R$ case, $\bar {\epsilon }_{11} \hat T_{0}$ could be greater than $\bar {\epsilon }_{01} \breve T_{0}$ in the near-wall region, as observed in figure 9.

If we go to the second order, the perturbation of the streamwise velocity and temperature of the 3-D travelling waves would also approach constants, which requires a Stokes layer to satisfy the no-slip condition. Balancing the unsteady and viscous terms, we obtain that the thickness of the Stokes layer is $y\sim R^{-1/2}$ for $\omega =O(1)$. Since they do not affect the leading-order balance, the Stokes-layer analysis is omitted in this paper; the detailed Stokes-layer solution can be found from Dong et al. (2020).

There is another issue that we would like to emphasise. From (4.28), we know that

(4.54)\begin{equation} \breve v_{01}\sim \hat C y\ln y, \quad\mbox{as }y\to 0, \end{equation}

where $\hat C$ is a constant. Under the wall-layer coordinate, this asymptotic behaviour is translated to

(4.55)\begin{equation} \breve v_{01}\sim \epsilon \ln \epsilon Y\hat C+\epsilon \hat CY \ln Y. \end{equation}

In the wall layer where $Y=O(1)$, the two terms on the right-hand side separate into two different scales. The wall-layer expansion (4.33) indicates that the above scale separation leads to two orders of solutions $(\breve U_1,\breve V_1,\breve W_1,\breve T_1)$ in the wall layer. The leading-order solution (4.41) matches the leading-order term in (4.55), where the $O(1)$ part of (4.41) induces an outflux to the main layer. To avoid lengthy mathematical argument, we do not show the second-order solution in the wall layer. Although ignoring the second-order solution of the wall layer will prevent us from constructing a consistent composite solution for the whole boundary layer, the leading-order solution is sufficient to derive a viscous correction to the main-layer solution, as will be demonstrated in the next subsection. Therefore, the numerical justification of the perturbation profiles, as will be displayed in figures 19 and 20, will be only focus on the comparison of the asymptotic predictions and the NPSE calculations in the main layer.

4.4. Improved asymptotic theory

Indeed, the second-order terms in the main-layer expansion (4.17) should be of $O(\epsilon \ln \epsilon )$, which is driven by the wall-layer outflux $\epsilon \ln \epsilon V_\infty$ and $\epsilon \ln \epsilon V_{1\infty }$. Since the value of $\epsilon \ln \epsilon$ is usually not small, neglecting it may lead to a quantitatively large error. As demonstrated by Dong et al. (2020), if the corrections from the lower boundary are taken into account, the accuracy of the main-layer solutions could be improved significantly. Thus, the improved boundary conditions for (4.20d–f) are derived,

(4.56)\begin{equation} \hat v_{01}(0)=\epsilon\ln\epsilon V_{\infty},\quad \hat v_{11}(0)=\hat v_{1-1}(0)=\epsilon \ln \epsilon V_{1\infty}, \end{equation}

where $V_\infty$ and $V_{1\infty }$ were defined in (4.43) and (4.53), respectively. Since $\hat A_3$ and $\hat A_2$ in $V_\infty$ and $V_{1\infty }$ are functions of $\hat p_{11}(0)$, the boundary condition (4.56) is homogeneous. Now the improved strategy is to solve the eigenvalue system (4.21) with boundary conditions (4.20a–c) and (4.56). The improved approach includes the impact of Reynolds number explicitly.

4.5. Discussion

From the above asymptotic analysis, we have described the skeleton of the fundamental resonance in hypersonic boundary layers by a triad resonance system appearing among the 2-D fundamental mode $(1,0)$, the 3-D travelling mode $(1,1)$ and the streak mode $(0,1)$. Here, the most distinguished feature is that the amplitude of the streamwise velocity component of the streak mode (streak component) is much greater than those of the transverse and lateral velocity components (roll components); therefore, the magnitude of the terms in the momentum equation governing the streak component is much greater than that governing the roll components. The mutual interaction of modes $(0,1)$ and $(1,1)$ could drive the formation of the roll components (see (4.19c,d)), but is too small to affect the leading-order streamwise momentum equation of the streak mode.

The conditions for the triad resonance include: (1) the dimensionless growth rates of the streak mode $(0,1)$ and the 3-D travelling modes ($1,\pm 1$) are of the same order as the dimensionless amplitude of the fundamental mode $\bar {\epsilon }_{10}$ (see (4.15a)); (2) the magnitude of the roll structure is smaller by a factor of $O(\sigma )$ than that of the streak structure (see (4.12)). Such scaling relations could not be prescribed by the SIA, which can also be regarded as a distinguished feature of the FR regime.

The value of the present asymptotic analysis is to reveal the intrinsic mechanism of the fundamental resonance from the dynamic point of view, as summarised in the schematic in figure 10. The nonlinear interaction of the fundamental mode and the streak mode seeds the growth of the 3-D travelling mode; the nonlinear interaction of the fundamental Mack mode $\hat \phi _{10}$ and the 3-D travelling modes $\hat \phi _{1\pm 1}$ drives the roll components of the streak mode, $\hat v_{01}$ and $\hat w_{01}$; the stronger amplification of the streak component of the streak mode is due to the linear lift-up mechanism.

Figure 10. Sketch of the principle of the FR.

The aforementioned analyses also indicate that the transverse asymptotic structures for different Fourier components are different. As shown in figure 11, the 2-D fundamental mode $\hat \phi _{10}$ shows a double-deck structure, namely, a main layer of $O(1)$ and a thinner Stokes layer of $O(R^{-1/2})$; the streak mode $\hat \phi _{01}$ shows an asymptotic structure with a main layer of $O(1)$ and a wall layer of $O(\epsilon )$; for the 3-D travelling mode $\hat {\phi }_{11}$, a main layer of $O(1)$, a wall layer of $O(\epsilon )$ and a thinner Stokes layer of $O(R^{-1/2})$ are observed. The wall-layer solutions of $\hat {\phi }_{01}$ and $\hat {\phi }_{1\pm 1}$ communicate with the main-layer solutions via outflux velocities, which are both $O(\epsilon \ln \epsilon )$. Inclusion of these effects leads to construction of the improved boundary conditions of the main-layer equations, which could increase the accuracy of the asymptotic predictions of both the growth rates and the perturbation profiles in the main layer.

Figure 11. Asymptotic structures of (a) $\hat {\phi }_{10}$, (b) $\hat {\phi }_{01}$ and (c) $\hat {\phi }_{1\pm 1}$, where $\epsilon \equiv (\bar {\epsilon }_{10} R)^{-1/3}$.

5. Verification of the asymptotic theory by NPSE calculations for moderate $R$ values

5.1. Parameters of the case studies

To verify the asymptotic theory in § 4, we choose a set of case studies with the same Mach number but different wall temperatures and Reynolds numbers, as listed in table 3. Each case is labelled by a two-digit character, the first and second of which distinguish the wall temperature and the Reynolds number, respectively. The frequency $\omega _0$, spanwise wavenumber $\beta _0$ and initial amplitude $\epsilon _{10}$ of the fundamental 2-D mode are the same as those in § 3, but the initial amplitudes of the oblique modes $\epsilon _{1\pm 1}$ are reduced to be $2.5\times 10^{-9}$, which enlarges the streamwise region of the fundamental resonance.

Table 3. Parameters for case studies in § 5.

5.2. Evolution of the fundamental modes

Figure 12(a,c) show the amplitude evolution of the fundamental mode $(1,0)$ and the MFD $(0,0)$ obtained by NPSE calculations for Cases A1–A5 and Cases B1–B5, respectively, and figure 12(b,d) display their zoom-in plots in the nonlinear phase. For all the Reynolds numbers, the mode $(1,0)$ becomes saturated at $x\approx 1000$ for Case A and $x\approx 800$ forCase B. Overall, the saturated amplitudes of both $(1,0)$ and $(0,0)$ components decrease with increase of $R$, except mode $(1,0)$ for Case A1. Since the MFD is mainly driven by the self-interaction of the fundamental mode, its order of magnitude is square of that of the latter.

Figure 12. Evolution of the amplitude of modes $(1,0)$ and $(0,0)$ obtained by NPSE: (a,b) results for Cases A1–A5; (c,d) results for Cases B1–B5. Panels (b,d) are the zoom-in plots of panels (a,c), respectively.

Figure 13 compares the profiles of the streamwise velocity and temperature of the base flow $(U_B,T_B)$ with those of the mean flow $(\bar U,\bar T)$. The difference between the two families of curves is quite limited.

Figure 13. Comparison of the base flow $\varPhi _{B}$ and the time- and spanwise-averaged mean flow $\bar \varPhi$ for Case A5 and Case B5 at $x=1000$: (a) streamwise velocity; (b) temperature.

In figure 14, we compare the wavenumber and growth rate of the Rayleigh solutions and NPSE calculations, where only Cases A5 and B5 are chosen for demonstration. Although there exist small discrepancies between the two families of curves, the overall trends agree. The error is attributed to the nonlinear, non-parallel and viscous effects.

Figure 14. Comparison of (a) the wavenumber $\alpha _r$ and (b) the growth rate $-\alpha _i$ of mode $(1,0)$ obtained by the NPSE calculations and the Rayleigh solutions for Case A5 and Case B5.

Figure 15 shows the comparison of the perturbation profiles obtained by Rayleigh and NPSE calculations. The largest error appears in the near-wall region and the critical layer, because the viscosity, neglected in the Rayleigh equation (4.5), is more important there. However, the overall agreement is quite satisfactory.

Figure 15. Comparison of the perturbation profiles of the fundamental mode $(1,0)$ obtained by the NPSE calculation and the Rayleigh solutions for (a,b) Case A5 and (c,d) Case B5 at $x=1000$.

5.3. The 3-D travelling mode and the streak mode

Figure 16(a) shows the evolution of the amplitudes of (1,1) and (0,1) obtained by NPSE calculations for Cases A1–A5. Being similar to figure 7, a mild increase for the 3-D travelling mode $(1,1)$ before $x \approx 720$ is observed, which is determined by its linear instability. After that, a drastic amplification occurs for both $(1,1)$ and $(0,1)$. In panel (b), we also plot their growth rates, which are nearly identical and increases monotonically in the interval of $x \in [900,1200]$, indicating the FR phenomenon. The increase of the growth rates is attributed to the increase of the amplitude of the fundamental mode $(1,0)$, as predicted by the asymptotic analysis (4.15a). Similar observations can be found in panels (c,d) for Cases B1–B5. However, the position where FR appears is promoted to $x \approx 650$, because the fundamental mode reaches a finite amplitude at an earlier position due to its higher growth rate.

Figure 16. Evolution of the (a,c) amplitude and (b,d) growth rate of $(1,1)$ and $(0,1)$ components: (a,b) Cases A1–A5; (c,d) Cases B1–B5.

Comparing with figure 14(b), we find that the FR appears before the mode $(1,0)$ reaches the nonlinear saturation phase. At the FR regions ($x\in [900,1200]]$ for Case A and $x\in [650,1000]$ for Case B), both the amplitude of the fundamental mode $\bar \epsilon _{01}$ and the growth rate of the secondary instability mode $\sigma$ increase with $x$, agreeing qualitatively with the scaling relation (4.15a).

To confirm quantitatively the asymptotic prediction of the growth rate of the secondary instability mode, figure 17 compares the NPSE results with the asymptotic predictions of both the original and improved versions. The growth rate of the NPSE calculation is normalised by the amplitude of the fundamental mode at each $x$,

(5.1)\begin{equation} \bar \sigma^{(m,n)}(x)= \frac{\sigma_{E}^{(m,n)}(x)}{ \dfrac{1}{2} \tilde{u}_{max}^{(1,0)}(x) }, \end{equation}

where $\sigma _{E}^{(m,n)}$ and ${u}_{max}^{(1,0)}$ are defined in (3.5) and (3.3), respectively. Here, the growth rates of $(0,1)$ and $(1,1)$ obtained by NPSE are both shown. In the FR regions, $x\in [900,1200]$ for Case A5 and $x\in [650,1000]$ for Case B5, the normalised growth rates $\bar \sigma ^{(0,1)}$ and $\bar \sigma ^{(1,1)}$ remain almost constant, confirming the scaling relation (4.15a). The growth rate of the original asymptotic theory $\bar \sigma ^{{Asmp}}$ is obtained by solving (4.21) with (4.19), shown by the red lines in the figure. Overall, $\bar \sigma ^{{Asmp}}$ does not change much with $x$ in the FR region, showing the same feature as the NPSE calculations, but its value is almost three times greater than $\bar \sigma ^{(0,1)}$ or $\bar \sigma ^{(1,1)}$. If the impact of the viscous wall layer is taken into account, as predicted by the improved asymptotic theory by solving (4.21) with (4.20a–c) and (4.56), then the agreement with the NPSE results is much better, as indicated by the blue lines in the figure.

Figure 17. Dependence on $x$ of the normalised growth rate obtained by the NPSE calculation, the asymptotic prediction and the improved asymptotic prediction for (a) Case A5 and (b) Case B5.

In figure 18, we show the impact of $R$ on the predictions of the growth rate $\bar \sigma$, where only the curves for $\bar \sigma ^{(0,1)}$ are plotted for the NPSE calculations. Again, the improved asymptotic predictions show a better agreement with the NPSE calculations than the original asymptotic prediction, and its relative error is only approximately $20\,\%$. The discrepancy is acceptable because in the asymptotic theory, the impact of the high-order viscosity, the non-parallelism and the higher-order Fourier components are excluded.

Figure 18. Dependence on $R$ of the normalised growth rate $\bar \sigma$ obtained by the NPSE calculation, the asymptotic prediction and the improved asymptotic prediction at $x=1000$: (a) Cases A1–A5; (b) Cases B1–B5.

Figures 19 and 20 compare the perturbation profiles of the streak mode (0,1) and the 3-D travelling mode $(1,1)$ obtained by the NPSE calculations and the main-layer asymptotic predictions (4.18) and (4.19) for Cases A and B. According to the system (4.19), $\hat {u}_{01}$, $\hat {v}_{01}$ and $\hat {T}_{01}$ are almost real, but $\hat {w}_{01}$ is almost pure imaginary, and therefore, we display only $\mathrm {Re}(\hat {u}_{01})$, $\mathrm {Re}{\hat {v}_{01}}$, $\mathrm {Im} (\hat {w}_{01})$ and $\mathrm {Re}(\hat {T}_{01})$. Note that the asymptotic predictions in these figures are only from the main-layer equations (4.18) and (4.19), and the solutions in the wall layer shown in § 4.3 are not included. In principle, if we construct a composite solution based on both the main-layer and wall-layer solutions, then the agreement with the NPSE calculations should be throughout the whole boundary layer. However, as mentioned before, because the main-layer solution undergoes a logarithmic singularity as the wall is approached, a consistent composite solution requires consideration of the wall-layer solutions up to the second order, which requires more complicated mathematical processes but yields the same accuracy on predicting the growth rate. Thus, we only probe the leading-order wall-layer solution, which is sufficient to construct the improved asymptotic approach as shown in § 4.4. Therefore, in figures 19 and 20, we only consider the comparison in the main layer ($\kern0.7pt y>2$). The original asymptotic results shown by the dashed lines are not sufficient to predict the NPSE calculations, even when $R$ is as large as $1 \times 10^{4}$. The discrepancy is greater in the region of $y<10$. However, when the wall-layer-induced correction is considered, the accuracy of the asymptotic predictions is remarkably improved, as shown by the symbols. In the region around the critical layer, the viscous effect needs to be taken into account to achieve a better agreement. A similar observation is obtained for the comparison of $|\hat {\phi }_{11}|$, as shown in panels (e–h).

Figure 19. Perturbation profiles of $(0,1)$ and $(1,1)$ at $x=1000$ for Case A. Solid lines, NPSE results; dashed lines, asymptotic predictions; symbols, improved asymptotic predictions. The profiles are normalised by the maximum of $\mathrm {Re}(\hat {v}_{01})$.

Figure 20. Perturbation profiles of $(0,1)$ and $(1,1)$ at $x=1000$ for Case B. Solid lines, NPSE results; dashed lines, asymptotic predictions; symbols, improved asymptotic predictions. The profiles are normalised by the maximum of $\mathrm {Re}(\hat {v}_{01})$.

In the main layer, the magnitude of $\hat {u}_{01}$ is much greater than those of $\hat {v}_{01}$ and $\hat {w}_{01}$, showing a longitudinal streak nature, which agrees with the scaling relations in (4.12). From the NPSE calculations, as $R$ increases, the peak of $\mathrm {Re} (\hat {u}_{01})$ increases monotonically, with its location moving towards a wall, agreeing with the improved asymptotic predictions. The main-layer solution of $\hat {v}_{01}$ and $\hat {w}_{01}$ are almost independent of $R$, which also agrees with the improved asymptotic predictions. For the 3-D travelling mode $\hat {\phi }_{11}$, the discrepancy between the improved solution and the asymptotic one is even smaller than that for $\hat {\phi }_{01}$, and both methods can provide a satisfied prediction in comparison with the NPSE calculations.

6. Verification of the asymptotic theory by the SIA approach for large $R$ values

In this section, we verify the asymptotic predictions for large $R$ values. Because the NPSE calculations for higher Reynolds numbers are not stable numerically, we verify our asymptotic predictions for higher $R$ values by the SIA approach in this section, which is confirmed to be accurate at moderate $R$ values shown in figures 8 and 9. The comparison is made based on a wavy base flow, which consists of a laminar Blasius solution and a fundamental mode $(1,0)$ with a finite amplitude. The latter is obtained by solving the O-S equation (2.10). Admittedly, this base flow is a bit artificial, but it is easy to demonstrate the fundamental resonance for sufficiently high Reynolds numbers. The parameters of the case studies are listed in table 4.

Table 4. Parameters for the FR calculations at higher $R$ values.

Figure 21 shows the perturbation profiles of the fundamental mode $(1,0)$ for Cases A and B with different $R$ values. The results of the Rayleigh equation (4.5) are also plotted by the dashed lines for comparison. Obviously, the planar Mack second mode $(1,0)$ shows a double-deck structure, a main layer in the major part of the boundary layer, where the O-S and Rayleigh solutions agree, and a Stokes layer in the near-wall region, where the streamwise velocity and temperature damp to satisfy the no-slip and isothermal conditions. An overshoot is observed at the edge of the Stokes layer for each curve, which can be predicted by the Stokes-layer solution. The thickness of the Stokes layer $\delta _S$ decreases as $R$ increases, and a scaling of $\delta _S \sim R^{-1/2}$ is evident. The eigenfunction $|\hat {T}_{10}|$ shows a sharp peak at the critical layer but $|\hat {u}_{10}|$ does not. This is because for a 2-D case, the streamwise velocity $|\hat {u}_{10}|$ shows only a logarithmic singularity, much weaker than the first-order singularity of $|\hat {T}_{10}|$. It is found that as $R$ increases, the O-S solutions approach the Rayleigh solution, except in the near-wall region, confirming the inviscid nature of the (1,0) component.

Figure 21. Perturbation profiles of the fundamental mode $(1,0)$ for the cases listed in table 4: (a,c) streamwise perturbation velocity; (b,d) perturbation temperature. (a,b) Cases A6–A10; (c,d) Cases B6–B10.

In figure 22, we compare the growth rate $\bar \sigma$ of the secondary mode obtained by the SIA approach, the asymptotic theory and the improved asymptotic theory. The amplitude of the fundamental mode $(1,0)$ is chosen to be $\bar {\epsilon }_{10}=0.015$. The growth rate $\bar \sigma$ obtained by the SIA increases with $R$ monotonically and approaches a constant as $R \to \infty$, indicating that the viscosity plays a stable role. In the large-$R$ limit, the relative errors between the SIA solutions and the asymptotic predictions are 18 % and 28 % for Cases A and B, respectively. Although there is still a gap between the SIA solution and the asymptotic prediction when $R$ is as high as $10^{8}$, the convergent trend is clearly seen. Including the wall-layer correction to the asymptotic theory, the improved asymptotic theory leads to an increase of the accuracy at moderate $R$ values. The monotonic increase of $\bar \sigma$ with $R$ for SIA is also reproduced by the improved asymptotic theory. Actually, in our asymptotic analysis, there exist two small parameters, namely, $R^{-1}$ and $\bar \epsilon _{10}$, and the $O(R^{-1/2})$ and $O(\bar \epsilon _{01}^2)$ terms are neglected in the governing equations. In figure 22, although we have probed the results for large $R$ values, the amplitude of mode (1,0) $\bar \epsilon _{01}$ is kept unchanged, which is the reason why there still exists a discrepancy even when $R=10^8$. To confirm this argument, we perform calculations of the SIA by changing the amplitude of the fundamental mode $\bar \epsilon _{01}$ with fixed Reynolds numbers, as shown in figure 23. For each Reynolds number, the normalised growth rate obtained by SIA increases with the decrease of $\bar \epsilon _{01}$ and approaches a constant when $\bar \epsilon _{01}$ is sufficiently small. The constant agrees well with the normalised growth rate predicted by the improved asymptotic theory.

Figure 22. Dependence on $R$ of the normalised growth rate $\bar \sigma$ obtained by the SIA approach, the asymptotic theory and the improved asymptotic theory: (a) Cases A6–A10; (b) Cases B6–B10.

Figure 23. Dependence on ${\bar {\epsilon }}_{10}$ of the normalised growth rate $\bar \sigma$ obtained by the SIA and the improved asymptotic prediction: (a) Cases A6 and A10; (b) Cases B6 and B10.

Figure 24 compares the perturbation profiles of the streak mode $(0,1)$ obtained by the SIA and the asymptotic predictions for different $R$ values for Cases A6–A10. As $R$ increases, the SIA solutions approach consistently the asymptotic predictions for all the perturbation quantities. Similar observations can be found for Cases B6–B10, as shown in figure 25. The agreement of both the growth rate and the perturbation profiles between the SIA and the improved asymptotic predictions in the large-$R$ limit confirms the accuracy of our asymptotic analysis.

Figure 24. Comparison of the perturbation profiles for streak mode $(0,1)$ obtained by SIA (for Cases A6–A10) and asymptotic theory. The eigenfunctions are normalised by the maximum of $\mathrm {Re}({\hat v}_{01})$.

Figure 25. Comparison of the perturbation profiles for streak mode $(0,1)$ obtained by SIA (for Cases B6–B10) and asymptotic theory. The eigenfunctions are normalised by the maximum of $\mathrm {Re}({\hat v}_{01})$.

7. Concluding remarks and discussion

In this paper, we focus on the fundamental resonance in hypersonic boundary layers, appearing when the planar Mack second mode dominates the nonlinear phase, and the infinitesimal oblique travelling waves and the streak mode amplify with a rather large rate. This regime is frequently observed in the laminar–turbulent transition in 2-D or axisymmetric hypersonic boundary layers since the planar Mack mode is always the most unstable linear instability mode in its laminar phase.

For a Mach 5.92 flat-plate hypersonic boundary layer with different wall temperatures and Reynolds numbers, we calculate the FR process using the NPSE approach, and show the amplitude evolution of representative Fourier components and their perturbation profiles. It is found that when the fundamental 2-D mode reaches a finite amplitude, a series of infinitesimal Fourier components with the same spanwise wavenumbers grow with the same rate, much greater than that of the fundamental mode, and the streak mode attains the greatest amplitude among these components. Such a phenomenon can be predicted quantitatively by the SIA based on a base flow consisting of the time- and spanwise-averaged mean flow and the fundamental mode. However, the SIA is not sufficient to reveal the underlying mechanism determining the energy transfer among different Fourier components and to explain the stronger amplification of the streak mode.

Therefore, a large-$R$ asymptotic theory is developed, based on the weakly nonlinear framework. The asymptotic analysis indicates that the FR is in principle a triad resonance system appearing among a dominant planar fundamental mode, an oblique travelling mode with the same frequency as the fundamental mode and a streak mode with the same spanwise wavenumber as the oblique mode. Remarkably, the amplitude of the streamwise velocity component (streak component) of the streak mode is much greater than those of the transverse and lateral velocity components (roll components), implying that these components may be driven by different mechanisms. The triad resonance system sketched in figure 10 shows that in the major part of the boundary layer, the interaction of the fundamental mode and the streak mode seeds the growth of the oblique mode, whereas the nonlinear interaction of the fundamental mode and the oblique mode drives the roll component of the streak mode, which further encourages a stronger amplification of the streamwise component of the streak mode due to the linear lift-up mechanism. The triad resonance system appears when: (1) the dimensionless growth rates of the streak mode and the oblique mode are of the same order of the dimensionless amplitude of the fundamental mode $\bar \epsilon _{10}$; and (2) the amplitude of the streak mode is $O(\bar \epsilon _{10}^{-1})$ greater than the oblique mode. These observations indicate that the present asymptotic theory is superior to the SIA by providing an in-depth understanding of theFR mechanism.

The asymptotic analysis also reveals the multi-layered structure of the perturbations in the FR regime as sketched in figure 11. The main-layer solutions of the streamwise velocity, spanwise velocity and temperature of both the streak mode and the oblique mode become singular as the wall is approached, and hence a viscous wall layer needs to be taken into account. It is found that the wall layer produces an outflux velocity of $O(\epsilon \ln \epsilon )$ to the main layer, inclusion of which leads to an improved asymptotic theory. Comparing with the NPSE calculations for moderate Reynolds numbers and the SIA for sufficiently large Reynolds numbers, it is found that the asymptotic theory can predict both the overall growth rate and the main-layer profiles of the infinitesimal perturbations, and the improved asymptotic theory could increase the accuracy of the growth-rate predictions remarkably. Additionally, for a moderate amplitude of the fundamental mode $\bar \epsilon _{01}$, the error of the improved asymptotic prediction does not vanish even when $R$ is sufficiently high; this is because the $O(\bar \epsilon _{01})$ terms are also neglected in the asymptotic analysis. Further decreasing $\bar \epsilon _{01}$ leads to a remarkable reduction of the error, confirming the accuracy of the asymptotic analysis in this paper.