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 (Reference Chuvakhov and Fedorov2016) 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. Reference Dong, Liu and Wu2020; Song, Zhao & Huang Reference Song, Zhao and Huang2020; Dong & Zhao Reference Dong and Zhao2021; Li & Dong Reference Li and Dong2021).

2.3.2. Nonlinear parabolised stability equations (NPSEs)

The NPSE approach (Bertolotti, Herbert & Spalart Reference Bertolotti, Herbert and Spalart1992; Chang & Malik Reference Chang and Malik1994) 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. (Reference Zhao, Zhang, Liu and Luo2016), and our code validation is provided in the appendix of Song et al. (Reference Song, Dong and Zhao2022).

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 (Reference Sivasubramanian and Fasel2015), Chen et al. (Reference Chen, Zhu and Lee2017) and Hader & Fasel (Reference Hader and Fasel2019), 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. Reference Chen, Zhu and Lee2017; Xu et al. Reference Xu, Liu, Mughal, Yu and Bai2020), and our code validation and discretisation method are provided by Song, Zhao & Dong (Reference Song, Zhao and Dong2023).

3.1. Case studies

For demonstration of the FR, we select a wind-tunnel condition of Maslov et al. (Reference Maslov, Shiplyuk, Sidorenko and Arnal2001), 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. (Reference Dong, Liu and Wu2020) and Dong & Zhao (Reference Dong and Zhao2021). 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 (Reference White2006, 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 Reference Mack1987). 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 Reference Dong, Zhang and Zhou2008). 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 (Reference White2006, p. 553).

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

(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 (Reference Sivasubramanian and Fasel2015), Chen et al. (Reference Chen, Zhu and Lee2017) and Hader & Fasel (Reference Hader and Fasel2019). 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. Reference Chen, Zhu and Lee2017; Hader & Fasel Reference Hader and Fasel2019). 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.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. (Reference Dong, Liu and Wu2020) and Dong & Zhao (Reference Dong and Zhao2021), 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. (Reference Dong, Liu and Wu2020), Dong & Zhao (Reference Dong and Zhao2021) and Zhao & Dong (Reference Zhao and Dong2022), 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 (Reference Schmid and Henningson2001, pp. 44–45); an example can be found from Dong et al. (Reference Dong, Liu and Wu2020). 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. (Reference Dong, Liu and Wu2020) 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 (Reference Schmid and Henningson2001, 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. (Reference Song, Dong and Zhao2022).

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 (Reference Malik1990) 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 (Reference Abramowitz and Stegun1964, 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. (Reference Dong, Liu and Wu2020).

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. (Reference Dong, Liu and Wu2020), 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.

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.