Interaction and breakdown induced by multiple optimal disturbances in hypersonic boundary layer

Abstract

The present paper finds that the coexistence of multiple primary instability waves may cause a non-trivial nonlinear interaction and breakdown process, which has not been reported before. In the considered Mach 6 flat-plate boundary layer, a global resolvent analysis reports three optimal disturbances (local maxima): a high-frequency planar wave, a low-frequency oblique wave and a stationary streak. For the dominant planar and oblique waves, a parabolised stability equation analysis identifies the initial non-modal transient growth and downstream modal growth. Initiated by these two optimal disturbances jointly, the complete linear and nonlinear instability processes until breakdown to turbulence are shown with direct numerical simulation. Owing to the transient growth, the oblique wave may be more significant than the planar wave in the breakdown. The oblique wave and scales of nonlinear interactions are pronounced in the outer layer, whose significance may not be comprehensively characterised by the wall pressure measurement. Fourier modes characterising the oblique-wave oblique breakdown, the planar-wave fundamental resonance, the planar-wave subharmonic resonance and the combination resonance related to a detuned mode are observed successively. The detuned mode seems to dominate the near-wall dynamics in the late nonlinear stage, characterised by $\varLambda$-like structures. Meanwhile, the existence of this detuned mode is independent of the initial amplitude ratio and the absolute amplitude of the oblique and planar waves. Weakly nonlinear stability analyses demonstrate that the detuned mode is mainly a consequence of the secondary instability under the combination of planar and oblique primary waves. Wave vector plots reveal the resonant state of multiple triads. Energy budget and amplitude-correction analyses provide a clear physical image of energy transfer.

1. Introduction

The boundary layer transition is of great significance in both flow physics and engineering applications. Many associated questions, such as what physical effect is dominant in the transition and how a flow breaks down to turbulence, are not fully understood. It is currently known that the receptivity, transient growth, eigenmode growth, parametric instability and resonance and bypass mechanisms can be involved in the transition process (Morkovin, Reshotko & Herbert 1994). Early knowledge of the boundary layer transition was partly provided by linear eigenmode analysis of Tollmien–Schlichting instability waves. However, in hypersonic states, a family of additional higher modes of acoustic nature were found to be important by linear stability theory (LST) with the parallel-flow assumption (Mack 1965, 1984). Particularly, the Mack second mode has received much attention due to its high growth rate if the Mach number exceeds approximately 4.

Despite the success of linear stability analysis, our understanding of the mechanisms of the boundary layer transition is still insufficient. Although the second mode appears to be preponderant in the linear instability stage, the upstream receptivity and transient growth as well as the downstream nonlinear stage are less considered jointly. In the flow past a blunt body, the initial amplitude of the first-mode waveband is potentially able to overtake that of the second mode if the receptivities across the bow shock wave (Fedorov et al. 2013), entropy layer (Goparaju & Gaitonde 2021) or roughness (Kuester & White 2015) are considered. Meanwhile, the integrated N-factor of the first mode can exceed that of the second mode (Masad 1993). Direct numerical simulation (DNS) on a Mach 15 flow over a parabola by Zhong (2001) showed that the maximum receptivity coefficient of the first mode was over twice that of the second mode in a wide frequency range of the forcing. In experimental research, fluctuations in a conventional wind tunnel are usually the most energetic in a relatively low-frequency range (Schneider 2001). A recent comparative study by Hildebrand et al. (2022) of various hypersonic quiet or noisy wind-tunnel environments showed that the power spectral density of background noise is the largest at a frequency f = O(1) kHz and exponentially decays at higher frequencies. These findings indicate that the first mode may not be weak in the vicinity of the leading edge. Therefore, due solely to the upstream effect, it is already difficult to determine beforehand whether the first or the second mode is dominant in various geometric configurations and free-stream conditions.

Considering the downstream nonlinear instability stage, the complete breakdown mechanisms and relative significance of the first and second modes in hypersonic states remain unclear. The complicated parametric resonance in a real environment is usually decomposed into individual regimes, including fundamental, subharmonic, oblique and combination resonance or breakdown. In a controlled hypersonic transition, the dominant primary wave is usually assumed or found to be the two-dimensional (2-D) second mode. Accordingly, numerous experimental (Chokani 1999; Kennedy et al. 2018a, 2022) and DNS (Sivasubramanian & Fasel 2015; Hader & Fasel 2019; Unnikrishnan & Gaitonde 2020) studies found that second-mode fundamental breakdown, with a frequency of hundreds of kilohertz, was dominant. In these observations, streamwise hot streaks with heat transfer overshoot and evident interactions between the second mode and its higher harmonics were common phenomena. Under the condition of a Boeing/AFOSR Mach 6 quiet tunnel, the second-mode fundamental breakdown was also significant in transitions induced by a wave packet (Sivasubramanian & Fasel 2014) or random forcing (Hader & Fasel 2018). As a lesser mechanism, second-mode subharmonic resonance was also detected (Bountin, Shiplyuk & Maslov 2008; Hader & Fasel 2019; Unnikrishnan & Gaitonde 2020). The second-mode oblique breakdown was numerically investigated by Franko & Lele (2013) and Hartman, Hader & Fasel (2021) but scarcely reported by experiments.

In other existing situations, however, low-frequency bands not belonging to the second-mode instability can be considerably energetic and participate in the nonlinear behaviour of multiple instability waves. These low frequencies could correspond to the first mode, Görtler vortices, wind-tunnel noise or the product of nonlinear interactions between high-frequency side bands. In addition to conventional ones such as the Görtler secondary instability (Ren & Fu 2015; Chen, Huang & Lee 2019) and first-mode oblique breakdown (Mayer, Wernz & Fasel 2007; Guo et al. 2022b), new breakdown scenarios involving these low-frequency components are possible. For example, a 20–30 kHz component was pronounced and even dominant on a flared cone in the Peking University Mach 6 quiet wind tunnel (Zhu et al. 2016; Xiong et al. 2020). To address this problem, Chen, Zhu & Lee (2017) employed the nonlinear parabolised stability equation (NPSE) to revisit the interaction between the low-frequency and the second mode, and indicated the possibility of combination resonance. However, further support for the corresponding experiments or full three-dimensional (3-D) DNS was not provided. Zhu et al. (2022) also experimentally demonstrated that a nearly adiabatic wall or heated wall enabled the first mode to be comparable to or overtake the second mode in frequency spectra. Under other wind-tunnel conditions, solid experimental evidence of interactions between the second mode and the low-frequency component was presented by Kimmel & Kendall (1991), Bountin et al. (2008), Munoz, Heitmann & Radespiel (2014) and Craig et al. (2019). Generally, however, the second mode was considerably stronger than the low-frequency component in the primary instability stage. The present paper serves to show another possibility that the low-frequency optimal disturbance can be as important as the optimal planar wave, which evolves into the second mode downstream. Furthermore, the present study aims to reveal the breakdown scenario and competition mechanism if the low-frequency component and the high-frequency second mode are equally significant and both result in primary instabilities.

In addition to the long-term discussions and debates on the first and second modes, non-modal instabilities have also attracted attention. Here, the terminologies ‘modal’ and ‘non-modal’ in the present paper denote the long-term characteristic of the eigenvalue problem for normal modes and the short-term characteristic of the initial value problem in local analyses, respectively (Schmid 2007; Kerswell 2018). With the appearance of numerous relevant studies in the 1990s, the physical image of transition had to be reconsidered. For instance, Trefethen et al. (1993) claimed that although eigenmodes of the linear Navier–Stokes (N–S) equation decayed monotonically, the non-modality of the operator led to an algebraic transient growth of the disturbances by factors of O(10⁵). Hanifi, Schmid & Henningson (1996) identified the ‘lift-up’ mechanism of non-modal growth in compressible boundary layers. Since then, remarkable accomplishments have been witnessed in linear and nonlinear non-modal theories (Schmid 2007; Kerswell 2018). One of the effective analytical tools is resolvent analysis, which seeks the maxima of the optimal gain in the frequency or wavenumber space. This approach has been intensively applied to identify the combination of the most responsive forcing and most amplified state in dynamic systems of wall-bounded shear flows (Ahmadi et al. 2019; Herrmann et al. 2021; Martini et al. 2021; Rigas, Sipp & Colonius 2021). Other recently popular analysis tools for optimal disturbances in spectral space are represented by the direct-adjoint PSE (Paredes et al. 2016) and the one-way N–S equation (Towne et al. 2022) approaches. These spatial marching methods require less memory usage than resolvent analysis based on matrix decomposition, which show advantages in cases with large numbers of grid nodes. The ‘one-way’ method further improves the robustness and accuracy for complicated problems compared with the PSE.

Based on the aforementioned novel methods, non-modal instabilities have been recognised to provide a crucial route to the high-speed boundary layer transition. Paredes, Choudhari & Li (2017, 2019) applied the NPSE to show that the optimal streak with proper spacing stabilised the first and second modes. Thus, interactions between modal and non-modal instabilities are of interest. Non-modal instability was also found to increase with higher Mach number (Tempelmann, Hanifi & Henningson 2012). In addition, non-modal instability is capable of evolving into the only dominator. Both NPSE and DNS (Paredes, Choudhari & Li 2020) as well as experimental (Kennedy et al. 2022) studies found that large nose bluntness could suppress the second mode and exhibit instabilities dominated by non-modal disturbances under flight conditions. The present paper is motivated to study the complete dynamics of multiple optimal disturbances in the linear, nonlinear and transitional regimes. The characterisation of early non-modal and subsequent modal instabilities of optimal disturbances is naturally included in a conventional framework of resolvent analysis and the PSE. The 3-D DNS is utilised to study the multiple instabilities in the late stage. The article is organised as follows. A description of flow conditions and numerical and theoretical tools is given in § 2. The computational set-up for the DNS, including the perturbation strategy of external forcings, is detailed in § 3. Results and discussions on the linear and nonlinear development of optimal disturbances are presented in § 4. Concluding remarks are provided in § 5.

2. Flow conditions and methodology

The flow condition of the Ludwieg-type short-duration wind tunnel Tranzit-M by Bountin et al. (2013) is considered here, whose hypersonic boundary layer stability and control has been widely studied. The free-stream conditions are given as follows: Mach number M _∞ = 6.0, stagnation temperature T ₀ = 354.08 K, static temperature T _∞ = 43.18 K and unit Reynolds number Re _1∞ = 1.05 × 10⁷ m⁻¹. The wall is assumed to be isothermal with a temperature of T_w = 293 K. The corresponding ratio of the wall temperature to the laminar-flow recovery value is T_w/T_rec = 0.954.

The flat-plate model for DNS has a length of L_x = 0.5 m with a sharp leading edge. A Cartesian coordinate system is constructed with the origin at the leading edge, the x direction along the direction of the free-stream velocity, the y direction perpendicular to the flat plate and the z direction satisfying the right-hand rule. The size of the effective computational domain for 3-D DNS is L_x × L_y × L_z = 0.5 m × 0.1 m × 0.01${\rm \pi}$ m, while additional sponge zones are placed near the outflow boundary to minimise the reflection of disturbances (Mani 2012). The spanwise length of the computational domain L_z is designed to be exactly four times the spanwise wavelength of the linearly optimal oblique wave for a convenient Fourier analysis, as discussed later.

2.1. Direct numerical simulation and base-flow solver

The 3-D compressible N–S equations in the Cartesian coordinate system are written in the following conservation form:

(2.1)\begin{equation}\frac{{\partial {\boldsymbol{Q}}}}{{\partial t}} + \frac{{\partial {\boldsymbol{F}}}}{{\partial x}} + \frac{{\partial {\boldsymbol{G}}}}{{\partial y}} + \frac{{\partial {\boldsymbol{H}}}}{{\partial z}} = \boldsymbol{\mathcal{B}}{\boldsymbol{f}^{\prime}},\end{equation}

where Q = (ρ, ρu, ρv, ρw, ρE)^T denotes the vector of conservative variables and F, G and H are the vectors of inviscid and viscous fluxes. Detailed expressions can be found in Anderson (1995). Here, ρ is the density; u, v and w are the velocity components in the x, y and z directions, respectively; and E is the total energy per unit mass. The right-hand side forcing vector $\boldsymbol{f}^{\prime} = (\,{f^{\prime}_1},{f^{\prime}_2},{f^{\prime}_3},{f^{\prime}_4},{f^{\prime}_5})$ in (2.1) is introduced to trigger the boundary layer transition in unsteady simulations. Specifically in the present paper, the operator $\boldsymbol{\mathcal{B}}$ constrains f′ to the wall-normal profile only at x ₀ = 0.04 m, which looks like a Dirac delta function in the streamwise direction. This forcing vector is assumed to possess a harmonic form with a small amplitude, whose details are determined by resolvent analysis later. The assumption of a calorically perfect Newtonian fluid is made with a constant specific heat ratio of γ = 1.4. Sutherland's law is adopted to calculate the dynamic viscosity μ, while the thermal conductivity κ is computed with a constant Prandtl number of 0.72. In addition, the vector Q is assumed to be decomposed into a 2-D steady base flow and a small disturbance as

(2.2)\begin{equation}{\boldsymbol{Q}}(x,y,z,t) = {{\boldsymbol{Q}}_b}(x,y) + {\boldsymbol{Q}^{\prime}}(x,y,z,t),\end{equation}

where the subscript ‘b’ and the prime denote base-flow and perturbation variables, respectively. The laminar base flow, satisfying the 2-D N–S equation, is the initial condition of the 3-D DNS and for stability analyses. Meanwhile, the steady base flow is confirmed to converge via the same numerical scheme and mesh as those of the 3-D DNS. Below, unless otherwise stated, the local base-flow and perturbed primitive variables are non-dimensionalised by the corresponding free-stream base-flow quantities, except that pressure p is by ${\rho _\infty }u_\infty ^2$. The reference length scale L_ref = 1 mm is used for non-dimensionalisation.

Both the full 3-D DNS and laminar base-flow simulation are performed using an in-house multi-block parallel finite-volume solver called PHAROS (Hao, Wang & Lee 2016; Hao & Wen 2020; Hao et al. 2021). This solver has successfully resolved the 3-D instability problem of a hypersonic boundary layer over a double cone via DNS, obtaining excellent agreement with experimental results (Hao et al. 2022). The inviscid fluxes are calculated using the modified Steger–Warming scheme (MacCormack 2014), which is extended to a higher order using the seventh-order weighted essentially non-oscillatory scheme (Jiang & Shu 1996). The viscous fluxes are discretised by the second-order central difference scheme. The three-stage third-order total variation diminishing Runge–Kutta method is employed for time marching. The boundary conditions are given as follows: free-stream conditions are imposed on the far-field boundary, extrapolation is used on the outflow boundary and isothermal, no-slip and no-penetration conditions are employed on the wall boundary. For 3-D DNS, a periodic condition is used on the spanwise boundary.

2.2. Resolvent analysis

Starting from the linearised N–S equations with a specified base flow, resolvent analysis aims at finding the most amplified disturbances in the desired parametric and spatial domains. Throughout the paper, the terminology ‘global resolvent analysis’ refers to the recently popular framework of resolvent analysis (Bugeat et al. 2019), which resolves the disturbance shape on the x–y plane globally instead of at a ‘local’ station. A more detailed description of the algorithm can be found in the work of Bugeat et al. (2019). Substituting (2.2) into (2.1) and then subtracting the base-flow equation yields the following form:

(2.3) \begin{equation} \frac{{\partial {\boldsymbol{Q}^{\prime}}}}{{\partial t}} + \frac{{\partial {\boldsymbol{F}^{\prime}}}}{{\partial x}} + \frac{{\partial {\boldsymbol{G}^{\prime}}}}{{\partial y}} + \frac{{\partial {\boldsymbol{H}^{\prime}}}}{{\partial z}} = {\boldsymbol{N}}^{\prime} + \boldsymbol{\mathcal{B}}{\boldsymbol{f}^{\prime}}, \end{equation}

where N′ is the nonlinear higher-order term. The complete form of (2.3) can be found in Paredes (2014). Dropping the higher-order term of (2.3) yields the linearised N–S equations as

(2.4)\begin{equation}\frac{{\partial {\boldsymbol{Q}^{\prime}}}}{{\partial t}} + \frac{{\partial {\boldsymbol{F}^{\prime}}}}{{\partial x}} + \frac{{\partial {\boldsymbol{G}^{\prime}}}}{{\partial y}} + \frac{{\partial {\boldsymbol{H}^{\prime}}}}{{\partial z}} = \boldsymbol{\mathcal{B}}{\boldsymbol{f}^{\prime}}.\end{equation}

Equation (2.4) can be reformulated as

(2.5)\begin{equation}\frac{{\partial {\boldsymbol{Q}^{\prime}}}}{{\partial t}} = \boldsymbol{\mathcal{A}}{\boldsymbol{Q}^{\prime}} + \boldsymbol{\mathcal{B}}{\boldsymbol{f}^{\prime}},\end{equation}

where $\boldsymbol{\mathcal{A}}$ is the linearised N–S operator related to the Jacobian matrices.

Bi-Fourier transforms can be performed in temporal and spanwise directions, which allows a solution Q′ to be written in the following form:

(2.6)\begin{equation}{\boldsymbol{Q}^{\prime}}(x,y,z,t) = \hat{{\boldsymbol{Q}}}(x,y)\,\textrm{exp(i}\beta z - \textrm{i}\omega t) + \textrm{c}\textrm{.c}\textrm{.}\end{equation}

Here, $\hat{{\boldsymbol{Q}}}$ is the complex eigenfunction, β is the spanwise wavenumber, ω is the real angular frequency and c.c. denotes the complex conjugate. Since Q′ is the linear response of the forcing, f′ can be written in a similar form:

(2.7)\begin{equation}{\boldsymbol{f}^{\prime}}(x,y,z,t) = \hat{{\boldsymbol{f}}}(x,y)\,\textrm{exp(i}\beta z - \textrm{i}\omega t) + \textrm{c}\textrm{.c}\textrm{.}\end{equation}

Substituting (2.6) and (2.7) into (2.5) gives

(2.8a,b)\begin{equation}\hat{{\boldsymbol{Q}}} = \boldsymbol{\mathcal{RB}}\hat{{\boldsymbol{f}}},\quad \boldsymbol{\mathcal{R}} = {( - \textrm{i}\omega \boldsymbol{\mathcal{I}} - \boldsymbol{\mathcal{A}})^{ - 1}},\end{equation}

which reflects the relationship between the forcing and its linear response (Bugeat et al. 2019). Here, $\boldsymbol{\mathcal{I}}$ is the identity operator.

In resolvent analysis, the purpose is to seek the forcing and response that maximise the energy amplification, i.e. the optimal gain σ defined by

(2.9)\begin{equation}{\sigma ^2}(\beta ,\omega ) = \mathop {\max }\limits_{\hat{{\boldsymbol{f}}}} \frac{{{{\|{\hat{{\boldsymbol{Q}}}} \|}_E}}}{{{{\|{\boldsymbol{\mathcal{B}}\hat{{\boldsymbol{f}}}} \|}_E}}}.\end{equation}

Here, Chu's energy (Chu 1965) is used to calculate the energy norm as

(2.10)\begin{equation}{\|{\hat{{\boldsymbol{Q}}}} \|_E} = \iint_\varOmega {{{\boldsymbol{Q}}^ \ast }\boldsymbol{\mathcal{M}}\hat{{\boldsymbol{Q}}}\,\textrm{d}\kern 0.06em x\,\textrm{d}y} ,\end{equation}

where the asterisk denotes the conjugate transpose, $\varOmega$ is the domain for the resolvent analysis and $\boldsymbol{\mathcal{M}}$ is the weight operator given by Bugeat et al. (2019). Furthermore, dimensionless Chu's energy containing primitive variables serves as an important indicator of the disturbance evolution in the present study. In detail, the form of (2.10) at a local station rather than an overall integral is given by

(2.11)\begin{gather}{E_{Chu}}(x) = \frac{1}{2}\int_0^\infty {\left[ {\bar{\rho }({{u^{\prime}}^2} + {{v^{\prime}}^2} + {{w^{\prime}}^2}) + \frac{{\bar{T}}}{{\gamma Ma_\infty^2\bar{\rho }}}{{\rho^{\prime}}^2} + \frac{{\bar{\rho }}}{{\gamma (\gamma - 1)Ma_\infty^2\bar{T}}}{{T^{\prime}}^2}} \right]\textrm{d}y} ,\end{gather}

where all the quantities are dimensionless, as stated in § 2.1. This definition consists of the fluctuation kinetic energy and a positive definite thermodynamic energy, which is commonly adopted in non-modal stability analyses. In 3-D DNS, the definition of E_Chu is replaced by the maximum of (2.11) in the spanwise direction.

As shown by Sipp & Marquet (2013), Bugeat et al. (2019) and Dwivedi et al. (2019), the optimisation problem (2.9) can be transformed into an eigenvalue problem as

(2.12)\begin{equation}{\boldsymbol{\mathcal{B}}^\ast }{\boldsymbol{\mathcal{M}}^{ - 1}}{\boldsymbol{\mathcal{R}}^ \ast }\boldsymbol{\mathcal{MRB}}\,\hat{{\boldsymbol{f}}} = {\sigma ^2}\hat{{\boldsymbol{f}}}.\end{equation}

The operator $\boldsymbol{\mathcal{A}}$ is discretised utilising the same methods as for (2.1) except that the inviscid fluxes are computed by the modified Steger–Warming scheme near discontinuities detected by a modified Ducros shock sensor (Hendrickson, Kartha & Candler 2018) and a central scheme in smooth regions. The settings of the boundary conditions are identical to those in DNS. The resulting discrete eigenvalue problem is solved by ARPACK software for a given β and ω of the regular mode (Sorensen et al. 1996). Specifically, $\boldsymbol{\mathcal{R}}$ and its conjugate transpose are computed using the SuperLU library (Li et al. 1999). The largest eigenvalue σ corresponds to the optimal gain, and the associated eigenfunction represents the optimal forcing. The optimal response can be obtained via (2.8). Prior grid convergence of resolvent analysis is conducted, and the converged results are further compared with the DNS and PSE results. More details about the resolvent analysis solver and the associated validation cases are described in Hao et al. (2023).

2.3. Parabolised stability equation analysis

Driven by the need to develop an efficient and accurate tool for stability analyses, Herbert & Bertolotti (1987), Herbert (1988) and Bertolotti (1991) proposed the fundamental idea of linear PSE (LPSE) and NPSE. Within the scope of convective instabilities, they assumed that the modal shape function $\hat{\boldsymbol{\psi} }$ of the disturbance is streamwise slowly varying, while the rapidly distorted part is absorbed into the exponential wavefunction. In detail, the disturbance is decomposed into the following form with a finite truncation of series:

(2.13)\begin{equation}{\boldsymbol{\psi} ^{\prime}}(x,y,z,t) = \sum\limits_{m ={-} M}^M {\sum\limits_{n ={-} N}^N {{{\hat{{\boldsymbol{\psi} }}}_{\!mn}}(x,y)\,\textrm{exp}\left( {\textrm{i}\int_{{x_{0,\textit{P}}}}^x {{\alpha_{mn}}\,\textrm{d}\kern 0.06em\tilde{x}} + \textrm{i}n{\beta_0}z - \textrm{i}m{\omega_0}t} \right)} ,}\end{equation}

where the vector ψ = (ρ, u, v, w, T)^T and the superscript ‘T’ represents the transpose, ${\hat{\boldsymbol{\psi}}_{\!mn}}$ and α_mn are the shape function and complex streamwise wavenumber of the Fourier mode (m, n), respectively, and β ₀ and ω ₀ are the specified fundamental spanwise wavenumber and angular frequency, respectively. Moreover, x _0,P is the distance from the inflow station of the PSE to the leading edge, and M and N are the modal numbers kept in the truncated Fourier series. Specially, mode (0, 0) is physically the mean flow distortion (MFD). Mode (0, n) with a non-zero n usually corresponds to a longitudinal vortex mode, which can evolve into a streak downstream with the streamwise vorticity gradually damped.

Substituting (2.13) into (2.3), dropping the forcing vector and neglecting the terms of $O(1/Re_0^2)$, where Re ₀ = ρ_∞u_∞l ₀/$\mu_{\infty}$ and l ₀ = ($\mu_{\infty}$x _0,P/ρ_∞u_∞)^1/2, finally gives rise to the NPSE as

(2.14)\begin{equation}({\boldsymbol{\mathcal{L}}_0} + {\boldsymbol{\mathcal{L}}_1}){\hat{\boldsymbol{\psi}}_{\!mn}} + {\boldsymbol{\mathcal{L}}_2}\frac{{\partial {{\hat{\boldsymbol{\psi}}}_{mn}}}}{{\partial x}} + \frac{{d{\alpha _{mn}}}}{{\partial x}}{\boldsymbol{\mathcal{L}}_3}{\hat{\boldsymbol{\psi}}_{\!mn}} = {\boldsymbol{\mathcal{F}}_{mn}}.\end{equation}

The base-flow-related operators ${\boldsymbol{\mathcal{L}}_0}$, ${\boldsymbol{\mathcal{L}}_1}$, ${\boldsymbol{\mathcal{L}}_2}$ and ${\boldsymbol{\mathcal{L}}_3}$ arise from the effects of locally parallel flow, non-parallel base flow, non-local shape function and non-local streamwise wavenumber, respectively, and ${\boldsymbol{\mathcal{F}}_{mn}}$ is the nonlinear term. Dropping ${\boldsymbol{\mathcal{F}}_{mn}}$ in (2.14) yields the LPSE, and then only maintaining the local operator ${\mathcal{L}_0}$ will generate an eigenvalue problem corresponding to LST. Detailed expressions of (2.14) are given by Paredes (2014). The boundary condition of the shape function is given by the common Dirichlet type as

(2.15)\begin{equation}{\hat{u}_{mn}},{\hat{v}_{mn}},{\hat{w}_{mn}},{\hat{T}_{mn}} = 0,\quad \textrm{at}\ y = 0\ \textrm{and}\ y \to \infty ,\end{equation}

except for the wall-normal velocity for the MFD at infinity:

(2.16)\begin{equation}\textrm{d}{\hat{v}_{00}}/\textrm{d}y = 0,\quad \textrm{at}\ y \to \infty ,\end{equation}

to include the growth effect of the boundary layer thickness arising from the MFD. To close the physical problem, an iterative scheme is used for the streamwise wavenumber of non-MFD modes during spatial marching as

(2.17)\begin{equation}{\alpha _{mn,new}} = {\alpha _{mn,old}} - \textrm{i}\frac{1}{{{E_{mn}}}}\int_0^\infty {\bar{\rho }\left( {\hat{u}_{mn}^ \ast \frac{{\partial {{\hat{u}}_{mn}}}}{{\partial x}} + \hat{v}_{mn}^ \ast \frac{{\partial {{\hat{v}}_{mn}}}}{{\partial x}} + \hat{w}_{mn}^ \ast \frac{{\partial {{\hat{w}}_{mn}}}}{{\partial x}}} \right)\textrm{d}y} ,\end{equation}

where

(2.18)\begin{equation}{E_{mn}} = \int_0^\infty {\bar{\rho }({{|{{{\hat{u}}_{mn}}} |}^2} + {{|{{{\hat{v}}_{mn}}} |}^2} + {{|{{{\hat{w}}_{mn}}} |}^2})\,\textrm{d}y} .\end{equation}

The streamwise wavenumber of the MFD is not updated, however, to improve the numerical robustness. This indicates that the growth of the MFD is entirely reflected in the shape function. The iteration procedure (2.17) and algebraic solution of (2.14) continue until the residual norm of α_mn is less than 10⁻⁷. Subsequently, the procedure is marched to the next streamwise station until the specified end boundary.

Regarding the numerical implementation, our in-house code CHASES is employed, which integrates the LST, LPSE, NPSE and sensitivity analyses. A series of validation cases for LST and LPSE have been carefully compared with both theoretical (Guo et al. 2020, 2021, 2022a) and DNS (Guo et al. 2022b; Cao et al. 2023) results. The validation of the NPSE code is shown in Appendix A. In terms of the numerical scheme for PSE, the second-order backward Euler scheme is used for streamwise marching, while the fourth-order central difference method is applied for the wall-normal difference. The Vigneron technique is adopted to eliminate the numerical instability resulting from residual ellipticity (Vigneron, Tannehill & Rakich 1978). In the LST analysis which also provides the initial condition for the PSE, the global algorithm via the Chebyshev pseudo-spectral method and a local algorithm using the fourth-order compact finite difference are employed for wall-normal discretisation (Malik 1990). The LST eigenvalue problem is solved by the Intel math kernel library.

The presented PSE results have been examined by grid convergence studies. The initial conditions of the PSE are mostly optimal responses obtained from resolvent analysis at x _0,P = 0.045 m, downstream of x ₀ = 0.04 m where the optimal forcings are localised in resolvent analysis and DNS. Additionally, LST profiles are also used for initialisation in § 4.1 for comparison. Those newly generated Fourier modes due to nonlinear interactions are assumed to be initially zero in § 4.2 or imposed by Fourier-transformed DNS profiles to improve the accuracy in § 4.3.

3. Case description and simulation strategy

Provided the base flow is prepared, the research procedure starts from a global resolvent analysis to identify all local maxima of optimal gain, which correspond to multiple optimal disturbances. Subsequently, the profiles of optimal disturbances or forcings serve to initiate the PSE and 3-D DNS to reveal the linear and nonlinear behaviour of disturbances.

The grid size for the main computational domain of the 3-D DNS is 4001 × 251 × 97, whose x–y mesh is entirely identical to that of the 2-D steady simulation. The orthogonal structured mesh is uniformly distributed in the streamwise and spanwise directions. In the wall-normal direction, the mesh is stretched with a minimum spacing of 2 × 10⁻⁵ m on the wall, and the transitional and turbulent boundary layers are ensured to contain 150 to 210 grid cells inside. The dimensionless grid spacings are Δx ⁺ ≈ 1.88, $\Delta y_{min}^ + \approx 0.30$ and Δz ⁺ ≈ 4.92, evaluated by the friction velocity and wall kinematic viscosity in the fully developed turbulent region x = 0.37 m via either 3D-DNS statistics or the skin friction by the van Driest II formula for C_f (Franko & Lele 2013; Guo et al. 2022b). In the linear and early nonlinear stages, the mesh contains at least 45 grid cells for each wavelength of the optimal planar wave (also the second mode downstream) in the x direction and 24 cells exactly for each wavelength of the optimal oblique wave (also the first mode downstream) in the z direction. Grid convergence of the amplitude evolution of the optimal disturbances has been confirmed in these instability stages.

Figure 1 displays the contours of the base-flow velocity and density fields. The dashed line, corresponding to x ₀ = 0.04 m, represents the location where the forcings are added in the resolvent analysis and 3-D DNS. This type of forcing appears to be a Dirac delta function looking in the streamwise direction and a constantly activated perturbation in the wall-normal, spanwise and temporal dimensions.

Figure 1. Contours of (a) streamwise velocity and (b) density of the laminar base flow non-dimensionalised by the free-stream value. Dashed line: location where the optimal forcings are imposed on the base flow.

In the DNS, the optimal forcings f′ are directly added to the right-hand side of the N–S equations only at x ₀ = 0.04 m, which resembles those of the resolvent analysis. The DNS forcings consist of those corresponding to the optimal planar wave ${\boldsymbol{f}^{\prime\!\!}_p}$, the optimal oblique wave ${\boldsymbol{f}^{\prime\!\!}_o}$ and an additional background noise term ${\boldsymbol{f}^{\prime\!\!}_n}$. The resolvent analysis of these optimal disturbances is detailed in § 4.1, and the resulting parametric set-up for the perturbation strategy in DNS is briefly introduced here.

The mathematical form of the forcing in DNS is given by

(3.1a)\begin{gather}{\boldsymbol{f}^{\prime}}({x_0},y,z,t) = {{\boldsymbol{f}^{\prime}\!\!}_p}({x_0},y,z,t) + {{\boldsymbol{f}^{\prime}\!\!}_o}({x_0},y,z,t) + {{\boldsymbol{f}^{\prime}\!\!}_n}({x_0},y,z,t),\end{gather}

(3.1b)\begin{gather}{{\boldsymbol{f}^{\prime}\!\!}_p}({x_0},y,z,t) = {\epsilon_1}\,{\hat{{\boldsymbol{f}}}_{\!p}}({x_0},y)\,\textrm{exp(i}{\beta _p}z - \textrm{i}{\omega _p}t) + \textrm{c}\textrm{.c}\textrm{.,}\end{gather}

(3.1c)\begin{gather}{{\boldsymbol{f}^{\prime}\!\!}_o}({x_0},y,z,t) = {\epsilon_2}[\,{\hat{{\boldsymbol{f}}}_{\!o}}({x_0},y)\,\textrm{exp(i}{\beta _o}z - \textrm{i}{\omega _o}t) + {\hat{{\boldsymbol{f}}}_{ \!- o}}({x_0},y){\rm exp} ( - \textrm{i}{\beta _o}z - \textrm{i}{\omega _o}t)] + \textrm{c}\textrm{.c}\textrm{.,}\end{gather}

(3.1d)\begin{gather}{{\boldsymbol{f}^{\prime}\!\!}_n}({x_0},y,z,t) = {\epsilon_3}[2r(y,z,t) - 1].\end{gather}

The complete description of the parameters in (3.1) is provided later, and a general introduction to the simulation cases is given first. Settings of DNS cases for comparative studies are listed in table 1. Case A1 is considered as the baseline one, and case A2 is introduced to include the background noise effect in real environments, i.e. the term ${\boldsymbol{f}^{\prime\!\!}_n}$ is non-zero in (3.1). The noise term ${\boldsymbol{f}^{\prime\!\!}_n}$ has an identical form to the forcing f′ as a right-hand side term of N–S equations. For the baseline case A1 and cases A3–A6, the amplitude ${\epsilon_3}$ of the term ${\boldsymbol{f}^{\prime\!\!}_n}$ is set to zero. In comparison, the effect of background noise is involved in case A2, where ${\epsilon_3} = 0.1$. In (3.1d), $r \in [0,1]$ is the white noise signal provided by a random number generator. The chosen value of ${\epsilon_3}$ makes the norm of ${\boldsymbol{f}^{\prime\!\!}_n}$ at least one order of magnitude larger than those of ${\boldsymbol{f}^{\prime\!\!}_p}$ and ${\boldsymbol{f}^{\prime\!\!}_o}$ simultaneously. This set-up allows adequate selection of additional potentially unstable frequencies and wavenumbers and tends to approach the flow state in a noisy environment. In DNS, the white noise term is added every 50 time steps for case A2 to allow its convection in the meshed domain (Cao et al. 2023). This interval corresponds to an angular frequency of 20ω_p, where ω_p is the angular frequency of the optimal planar wave. We have also simulated another DNS case forced only by the same noise term. It is found that no transition to turbulence occurs upstream of x = 0.5 m. Therefore, the considered white noise intensity should be insufficient to trigger transition in the unstable region of the eigenmodes. The white noise only serves to add initial randomness to the laminar flow region of case A2.

Table 1. Computational set-up for DNS cases.

The DNS cases A1 and A2 for the main discussion in the following contents impose a specific initial amplitude of the forcings. Without loss of generality, additional cases A3–A6 are simulated to examine the effects of the amplitude ratio and absolute amplitude of the forcings. As shown in table 1, cases A3–A5 are used to clarify the effect of the amplitude ratio by adjusting the modulus of the amplitude parameters ${\epsilon_1}$ and ${\epsilon_2}$ in (3.1), which are specified later. The responses of the mass-flow-rate fluctuations slightly downstream at x = 0.045 m are as follows: close but unequal in case A3, where ${(\rho u)^{\prime}_{max,p}} = 2{(\rho u)^{\prime}_{max,o}} = 0.01$; planar wave dominates in case A4, where ${(\rho u)^{\prime}_{max,p}} = 10{(\rho u)^{\prime}_{max,o}} = 0.01$; and oblique wave dominates in case A5, where ${(\rho u)^{\prime}_{max,o}}\, = \,10{(\rho u)^{\prime}_{max,p}} = 0.01$. Another case with a weaker perturbation intensity is simulated to show the impact of the absolute amplitude, where ${(\rho u)^{\prime}_{max,p}} = {(\rho u)^{\prime}_{max,o}} = 0.002$ at x = 0.045 m. Other settings remain identical to those of case A1.

For further details, first, let the fundamental angular frequency ω ₀ and spanwise wavenumber β ₀ be defined by ω ₀L_ref/u _∞ = 0.1 and β ₀L_ref = 0.8, respectively. In (3.1), the following resolvent analysis gives ω_p = 10ω ₀ and β_p = 0 for the optimal planar wave, and ω_o = 3ω ₀ and β_o = β ₀ for the optimal oblique wave. Interestingly, ω_p = 10ω ₀, with dimensional frequency of f ≈ 125.8 kHz, is nearly the second-mode frequency corresponding to the maximum wall pressure fluctuation under the considered experiment condition by Bountin et al. (2013). Thus, the optimal frequency should be representative. For the baseline case A1 in table 1, the complex amplitudes are specified by ${\epsilon_1} = (1.554 \times {10^{ - 3}},\,2.1835 \times {10^{ - 7}})$ for the planar wave and ${\epsilon_2} = ( - 7.6119 \times {10^{ - 4}}, - 6.4655 \times {10^{ - 5}})$ for the oblique wave. The modulus of the input forcing amplitude ${\epsilon_1}$ or ${\epsilon_2}$ determines the output amplitude of optimal disturbances linearly, while the phase of the complex amplitude results in the initial phase of output disturbances. These amplitudes are determined beforehand by resolvent analysis and PSE, which would result in the following modal properties. For the optimal planar wave, the initial amplitude forces the maximum dimensionless mass flow rate fluctuation at x = 0.045 m to be ${(\rho u)^{\prime}_{max,p}} = 1 \times {10^{ - 2}}$ with a zero phase angle. In other words, the complex amplitude ${\epsilon_1}$ can be obtained by a linear relation as ${\epsilon _1} /1 \times {10^{ - 2}} = {\epsilon_{1,try}}/{(\widehat {\rho u})_{max,p,try}}$, where the subscript ‘try’ indicates a test run of resolvent analysis for linear amplitude scaling. Meanwhile, for the optimal oblique wave, the amplitude ratio to that of the planar wave is scaled based on the resolvent analysis with an equivalent 1:1 initial forcing norm (see § 2.2), which yields ${(\rho u)^{\prime}_{max,o}} \approx 1.01 \times {10^{ - 2}}$ downstream at x = 0.045 m. As mentioned in § 1, it is reasonable to impose an initial amplitude on the first mode no smaller than that of the second mode. Furthermore, for the oblique wave propagated with a negative wave angle, the relationship between its forcing shape $\hat{f}_{-o}$ and $\hat{f}_{o}$ obeys the symmetry condition described in Park & Park (2013). The profiles of the forcing shape and phase angle are shown in figure 2, where the phase angle is in the range of (−180°, 180°]. Part of the non-zero forcing is located outside the boundary layer, which is consistent with the usual features of non-modal disturbances (Bugeat et al. 2019). For the planar wave, as shown in figure 2(c), the energetic forcings are mostly above the relative sonic line for Mack second mode. This observation may indicate a difference in properties between non-modal and modal disturbances.

Figure 2. Profiles of the five components of the optimal forcing that are imposed at x = 40 mm for DNS: (a) modulus and (b) phase angle corresponding to the oblique wave, and (c) modulus and (d) phase angle corresponding to the planar wave. The plotted forcing $\hat{f}_{1}$ is non-dimensionalised by ρ _∞u _∞/L_ref, $\hat{f}_{2}$, $\hat{f}_{3}$ and ${\hat{f}_4}$ by ${\rho _\infty }u_\infty ^2/{L_{ref}}$, and $\hat{f}_{5}$ by ${\rho _\infty }u_\infty ^3/{L_{ref}}$. Dashed line in (a): location of the boundary layer thickness δ _0.99. Dashed-dotted line in (c): location of the relative sonic line of the second mode (angular frequency equals ω_p), which satisfies U = c_r − a. Here, U, a and c_r represent the base streamwise velocity, base sound speed and phase velocity of the Mack second mode.

4. Results

4.1. Linearly optimal disturbances

Global resolvent analysis is performed over a wide range of spanwise wavenumbers and angular frequencies to determine the most amplified disturbances. The computational domain is extended from x ₀ = 0.04 m to x_e = 0.2 m, which is almost the pre-transitional region in the 3-D DNS. Figure 3 displays the contour of the optimal gain in the spectral domain. Generally, three types of optimal disturbances are identified in this free-stream condition, namely the optimal planar wave (10ω ₀, 0β ₀) appearing on the line β = 0, the optimal oblique wave (3ω ₀, β ₀) and the streak (10⁻⁵ω ₀, 1.5β ₀). Note that further lowering the frequency does not change the optimal gain of the streak. Among the three optimal disturbances, the planar and oblique waves possess the largest optimal gain, which below are called the (Fourier) modes (10, 0) and (3, 1), respectively. Figure 4 compares the modal shape of the streamwise velocity between the LPSE in figure 4(a,c,e) and the resolvent analysis in figure 4(b,d,f). In addition, the effect of adding the streak is compared in figures 4(g) and 4(h) by NPSE. The initial amplitudes of optimal disturbances are based on the resolvent analysis with a 1:1 forcing amplitude ratio, as stated in §§ 2 and 3. The agreement between LPSE and the resolvent analysis is excellent. The incident long-wavelength and upright short-wavelength patterns of the oblique and planar waves are observed. Based on figures 4(g) and 4(h), the streak does not essentially alter the combined disturbance shape owing to the weak gain in figure 3. Therefore, it is sensible to exclude the linearly optimal streak in 3-D DNS and consider the optimal oblique and planar waves only. Furthermore, it is observed that the coexistence of multiple optimal disturbances gives rise to an outer oblique-wave-dominant layer and an inner planar-wave-dominant layer.

Figure 3. Contours of optimal gain in the parameter space of the spanwise wavenumber and the angular frequency. Three plotted circles: local maxima of optimal gain corresponding to three types of optimal disturbances. Here, ω ₀L_ref/u _∞ = 0.1 and β ₀L_ref = 0.8.

Figure 4. Streamwise velocity contours of optimal disturbances u′, normalised by each maximum, by LPSE for (a) oblique wave, (c) planar wave and (e) streak, by resolvent analysis for (b) oblique wave, (d) planar wave and (f) streak, and by NPSE with the resolvent-based initial amplitude ratio for (g) oblique wave + planar wave and (h) oblique wave + planar wave + streak. Solid line: location of the base-flow boundary layer thickness δ _0.99. The maximum u′ for normalisation is 0.061 (a,b), 0.039 (c,d), 0.012 (e,f), 0.16 (g) and 0.058 (h).

Figure 5 shows a quantitative evaluation of the evolution of the optimal disturbances. The N-factor based on the Chu's energy is defined by

(4.1)\begin{equation}{N_s} = 0.5\,\textrm{ln}({E_{Chu}}/{E_{Chu,0}}),\end{equation}

where E_Chu is defined by (2.11) and E_{Chu ,0} denotes the value at x = 0.04 m. Thus, the N-factor characterises the exponential amplification of Chu's energy, and the streamwise derivative dN_s/${\rm d}\kern 0.06em x$ indicates the corresponding growth rate. A satisfactory agreement is achieved between LPSE and the resolvent analysis. In the meantime, LPSE initialised by LST profiles provides the evolution of the first and second modes denoted by open symbols. In figure 5(b), the vertical dashed lines mark the locations starting from which the difference in dN_s/${\rm d}\kern 0.06em x$ between the closed and open symbols is within 1%. The two locations, at x ≈ 0.136 m for the planar wave and x ≈ 0.164 m for the oblique wave, categorise the regions into upstream non-modal instability and downstream purely modal instability. Compared with the planar wave, the oblique wave possesses a stronger non-modal transient growth and weaker modal growth. This tendency is consistent with existing knowledge that the second mode is mostly dominant in the linear instability stage of hypersonic transitions. However, due to the prominent transient growth, the oblique wave has only a slightly lower N-factor than the planar wave at x_e = 0.2 m. This result shows the possibility that the oblique first mode is as important as the second mode even in the linear instability stage. Without the transient growth effect, the open symbols illustrate that the N-factor of the planar wave exceeds that of the oblique wave by approximately 0.8 at x_e, leading to a five-times-larger ratio of Chu's energy. Furthermore, E_Chu of the oblique wave is 11 times that of the streak at x_e, indicating that the streak has a negligible impact.

Figure 5. Streamwise development of the (a) N-factor and (b) growth rate of Chu's energy dN_s/${\rm d}\kern 0.06em x$. Filled symbols, LPSE results initialised by resolvent-analysis profiles; open symbols, LPSE results initialised by LST profiles. Vertical dashed lines in (b): locations downstream of which the difference in dN_s/${\rm d}\kern 0.06em x$ between filled and open symbols is less than 1%, which can be viewed as the starting locations of purely modal instabilities for the planar and oblique waves.

4.2. The DNS results: instability waves and breakdown

Having examined the linear dynamics of optimal disturbances, we concentrate on the behaviour of the oblique and planar waves in the nonlinear regime. Firstly, the baseline case A1 and the case A2 with background noise are mainly discussed. Figure 6 shows the instantaneous and disturbed flow fields on the x–y plane. In addition to the outer oblique wave and inner planar wave, new streamwise flow scales characterising nonlinear sum or difference interactions appear upstream of x_e = 0.2 m. The newly produced streamwise scales and downstream breakdown seem to occur earlier in the outer boundary layer than in the inner layer, as illustrated in figure 6(a). This observation may indicate the importance of oblique waves in the breakdown. These multiple streamwise scales are also displayed in figure 6(b) and a Rayleigh-scatter-like visualisation (Zhu et al. 2016) in figure 6(e). Furthermore, the ‘rope-like’ structure of the density fluctuation, commonly observed for the second mode in experiments (Laurence, Wagner & Hannemann 2016; Kennedy et al. 2018b) and DNS (Pruett & Chang 1995; Unnikrishnan & Gaitonde 2020), is captured near the generalised inflection point in figure 6(d) and quantitatively shown in figure 6(f). Here, the generalised inflection point is the location where (∂/∂y)(ρ_b∂u_b/∂y) = 0. In figure 6(f), a wavelength of around 5 mm is identified as the second-mode wavelength, since it is approximately twice the local boundary layer thickness δ _0.99 = 2.25 mm. However, another evident long wavelength of 18 mm is observed, which should not belong to the second mode. This differs from the previous conclusion in the literature that the ‘rope-like’ structure is attributed to the second mode (Laurence et al. 2016). To determine the wavelength more accurately, narrow bandpass filtering is applied to the time series of ρ′, where the examined centre frequencies are ω_o = 3ω ₀ and ω_p = 10ω ₀. The corresponding streamwise wavelengths of the filtered data at x = 0.1 m are around 18.4 and 5.8 mm for the two frequencies, respectively. In § 4.3, figure 14(b) further demonstrates that the evident density fluctuation ρ′ near the rope-like structure is jointly contributed to by the oblique and planar waves. The different contributors to the rope-like structure in this paper and in the literature may arise from the differences in the geometry and flow environment. The present results here serve to show a possibility of the significance of low-frequency waves.

Figure 6. Flow fields on the plane z = L_z/2 for case A1 about (a) instantaneous u, full view, (b) instantaneous u, (c) disturbance u′, (d) disturbance ρ′ and (e) Rayleigh-scatter-like visualisation with a cut-off temperature T = 1.2, and (f) the disturbance ρ′ along the generalised inflection point in (d), with all non-dimensionalised by each free-stream base-flow value. Arrow in (a): the location where optimal forcings are imposed. Solid line in (d): location of the generalised inflection point. Notations in (b): the streamwise wavelengths of the oblique wave and of those generated by sum or difference nonlinear interactions.

The boundary layer also presents very different breakdown scenarios at various wall-normal heights, as shown in figure 7. The near-wall structures appear to be positive velocity streaks in figure 7(a). The inner layer consists of elongated structures with large velocity fluctuations first as an aligned pattern upstream and then as a staggered pattern downstream in figure 7(b). The outer layer mainly contains the staggered $\varLambda$-vortices with small velocity fluctuations in figure 7(c). Breakdown is found to have occurred at x > 0.35 m for all three slices. These observations may indicate the existence of fundamental, subharmonic and oblique breakdown in different streamwise or wall-normal locations. As a visualisation of vortices, figure 8 shows the Q-criterion iso-surface. The early spanwise distortion of the vortices and the formation of $\varLambda$-vortices and hairpin vortices are displayed.

Figure 7. Instantaneous dimensionless velocity fluctuation field u′ (compared with the base flow) on the x–z planes for case A1 at (a) y = 0.1 mm ( y/δ _0.99 ≈ 0.016 at x = 300 mm), (b) y = 1.3 mm ( y/δ _0.99 ≈ 0.21 at x = 300 mm) and (c) y = 3 mm ( y/δ _0.99 ≈ 0.49 at x = 300 mm).

Figure 8. The Q-criterion iso-surface (L_ref/u _∞)²Q = 10⁷ in the range 0.04 m < x < 0.3 m for case A1, coloured by the dimensionless pressure $p/{\rho _\infty }u_\infty ^2$.

To obtain a closer observation of the near-wall dynamics, figure 9 shows the instantaneous skin friction coefficient C_f for both cases A1 and A2. The streamwise travelling wave, oblique wave and large-scale $\varLambda$-vortex are seen sequentially. According to figure 9(b), the spanwise scale is mainly dominated by half the wavelength of the optimal oblique wave, i.e. λ_{z ,o}/2, which corresponds to the wavenumber 2β ₀. This is mainly contributed by the streak (0, 2) characterising the oblique breakdown, as shown by the Fourier analysis next. In the streamwise direction, the planar-wave wavelength λ_{x ,p} is first dominant and the oblique-wave wavelength λ_{x ,o} is then highly pronounced. This can be understood by a subsequent modal analysis in which the planar wave saturates earlier than the oblique wave. What appears to be new is that staggered $\varLambda$-vortices are formed upstream of the eventual breakdown with a streamwise wavenumber that is not an integer multiple either of the oblique one or of the planar one. Hence, these $\varLambda$-vortices are inferred to be products of wave–wave interactions. Furthermore, the wavelength of the $\varLambda$-vortex is approximately five times the planar wavelength, corresponding to an angular frequency of 2ω ₀ provided that phase locking occurs. The physical mechanisms behind this are focused on in the following sections. Finally, by comparing figures 9(b) and 9(c), the white noise effect is found only to affect local meticulous structures (e.g., in circled regions of figure 9(c)) and the group velocity of disturbances at approximately x > 0.22 m. This finding suggests that the background noise does not have an evident impact on the interactions of the optimal disturbances and their higher harmonics.

Figure 9. Instantaneous skin friction coefficient C_f of (a) full view for case A1, (b) zoomed-in view for case A1 (without the background noise effect) and (c) zoomed-in view for case A2 (with the background noise effect). The symbols λ_z,o and λ_x,o denote the spanwise and streamwise wavelengths of the oblique wave (3, 1), respectively, and λ_x,p and λ_x,(2,2) represent the streamwise wavelengths of the planar wave (10, 0) and mode (2, 2), respectively.

Statistical results are shown to identify the pre-transitional, transitional and fully developed turbulent regions. The results are ensured to be independent from the selection of temporal observation windows. Figure 10(a) shows that transition begins at x ≈ 0.13 m with a minimum C_f and ends at x ≈ 0.37 m with a C_f value collapsing onto the turbulent correlation. Moreover, a small deviation between cases A1 and A2 in C_f or the boundary layer thickness exists at approximately x > 0.26 m. This finding shows that the background noise effect is not very significant to the mean flow evolution. Figure 10(b) shows the streamwise velocity spectrum at x = 480 mm. The agreement with the law of ω ^−5/3 in the inertial subrange and the law of ω ⁻⁷ in the dissipation scales suggests the establishment of a fully developed turbulent state at this station.

Figure 10. Quantitative results of (a) mean skin friction coefficient and boundary layer thickness for cases A1 (thick line) and A2 (thin line) and the van Driest II formula for C_f (dashed line), and (b) streamwise velocity spectrum ${E_{11}} = \langle u^{\prime}u^{\prime}\rangle /u_\infty ^2$ at x = 480 mm for case A1, where $\langle \cdot \rangle$ denotes the ensemble average.

The bi-Fourier transform into the t–z space gives information on a series of Fourier modes (m, n). Convergence of spectra in the pre-transitional and transitional regions has been achieved based on the following settings: Hann window function is used with 50% overlapping, and the lowest angular frequency is 0.05ω_p in each window. Firstly, for cases A1 and A2, figures 11(a) and 11(b) show the evolution of the maximum mass flow rate fluctuation and Chu's energy, respectively. The rapidly growing modes at x < 0.1 m include the oblique wave (3, 1) and the resulting streak (0, 2), the planar wave (10, 0), MFD (0, 0) and mode (10, 2) associated with the fundamental resonance. For these significant modes, the background noise only has a limited effect in the late stage of nonlinear instability and in the turbulent region. It is concluded that the background noise can neither select new significant modes with different frequencies or wavenumbers nor manipulate the main characteristics of the existing modes. In addition, modes (2, 0) and (2, 2) corresponding to the large-scale $\varLambda$-vortices in figure 9 and mode (5, 1) associated with the subharmonic resonance with the optimal planar wave start to be amplified rapidly in the vicinity of x = 0.1 m. In terms of the $\varLambda$-vortices, some neighbouring modes such as (1, 0) and (1.5, 0), with approaching streamwise scales to modes (2, 0) and (2, 2), are shown not to be pronounced in figure 11. Therefore, there should be a frequency selection process, which is discussed next. For brevity, evolution of some other significant modes is shown in Appendix B. The main interaction occurring in the oblique breakdown related to mode (3, 1), the fundamental resonance related to mode (10, 0) and the subharmonic resonance related to mode (10, 0) can be expressed by

(4.2)\begin{gather}(3,1) - (3, - 1) \to (0,2),\end{gather}

(4.3)\begin{gather}(10,\textrm{0}) + (0,\textrm{2}) \to (10,2),\end{gather}

(4.4)\begin{gather}(10,\textrm{0}) - (5,\textrm{1}) \to (5, - \textrm{1}),\end{gather}

respectively. The corresponding reverse interactions among the triads are also included, such as (0, 2) + (3, −1) → (3, 1) for the oblique breakdown.

Figure 11. Streamwise development of disturbances on (a) maximum mass flow rate ${(\rho u)^{\prime}_{max}}$ for cases A1 (solid line) and A2 (dashed line), (b) Chu's energy E_Chu for cases A1 (solid line) and A2 (dashed line), (c) Chu's energy E_Chu for case A4 (initial forcing ratio of planar to oblique waves equals 10:1), (d) Chu's energy E_Chu for case A5 (initial forcing ratio of planar to oblique waves equals 1:10) and (e) contour for logarithmic fast Fourier transform of wall pressure fluctuation for case A1, ln(|p′|).

The successive appearance of highly amplified modes (3, 1), (0, 2), (10, 0), (10, 2) and detuned (2, 2) or (2, 0) in figures 11(a) and 11(b) indicates the occurrence of oblique and fundamental breakdowns in conjunction with an unconfirmed combination resonance. With regard to the modal amplitude, the streak mode (0, 2) is dominant in the early region x < 0.2 m. Comparing the downstream first and Mack second modes, those of the oblique wave (3, 1) are stronger than and saturate later than those of the planar wave (10, 0). These results are consistent with the observations of flow structures in figures 6–9. However, the wall pressure fluctuation shows a considerable magnitude for the second-mode angular frequency 10ω ₀ at x < 0.2 m, as shown in figure 11(e). It can thus be inferred that the wall pressure measurement in experiments may not completely represent the possible modal competition in the outer layer of the boundary layer away from the wall, particularly for non-acoustic disturbances. The integral E_Chu, including the region away from the wall and containing kinetic and thermodynamic energy simultaneously, may characterise the modal growth and structure instabilities in figures 6–9 more comprehensively. In addition to the sum angular frequency 13ω ₀ and difference angular frequency 7ω ₀ in figure 11(e), the interaction between (3, 1) and (10, 0) generates a series of angular frequencies that are only integer multiples of the fundamental one ω ₀. These are caused by the fact that the two identified optimal frequencies 10ω ₀ and 3ω ₀ under this freestream condition are coprime. Physically, the spectral broadening during the transition could be enhanced by these discrete modes.

In terms of a comparison with other well-known breakdown scenarios, further information is shown by figure 11(c) for the planar-wave-dominant case A4 and by figure 11(d) for the oblique-wave-dominant case A5. Several interesting aspects are discussed as follows.

1. (i) For case A4 where the initial forcing ratio of planar to oblique waves equals 10:1, the planar wave (10, 0) saturates early in the vicinity of x = 0.2 m before the MFD is strong enough. Meanwhile, modes (10, n) are found not to be amplified enough, including those unshown with n ≠ 2. This indicates that fundamental resonance does not trigger the eventual breakdown due to a relatively weak second-mode wave in the considered flat-plate boundary layer, which may differ from those scenarios for a cone by Fasel's group. Particularly, the flared-cone flow adding the background noise by Hader & Fasel (2018) supported a constantly growing second mode with a slowly varying boundary layer thickness, which resulted in the spectral peak of the second mode at 300 kHz and its harmonics. Under our conditions, the optimal planar-wave frequency f ≈ 125.8 kHz, which is nearly the frequency of peak pressure fluctuations in the experiment by Bountin et al. (2013), does not possess a long unstable streamwise region. By contrast, the oblique wave of case A4 undergoes a much longer distance of the instability stage in this state. The amplified oblique wave exceeds the planar wave finally even with a lower initial amplitude. The downstream scenario in case A4 is more likely to be a joint type of oblique breakdown, combination resonance and second-mode subharmonic breakdown. The significance of low-frequency oblique waves in the late stage is somewhat consistent with the conclusions under the flow and geometric conditions at Peking University (Zhu et al. 2016; Chen et al. 2017).

2. (ii) The similarities between cases A1 and A5 upstream of x = 0.2 m suggest that the oblique breakdown is quite important in the early stage. Further downstream, with the existence of mode (10, 0), the possible detuned mode (2, 2) can be promoted. Among the comparative cases, case A1 reports the most prominent amplification of mode (2, 2). With regard to the signature of oblique breakdown, the streak (0, 2) of case A1 becomes weaker than that of case A5 in the saturation stage downstream of x = 0.2 m. Consequently, oblique breakdown plays a crucial role mainly in the early stage.

3. (iii) By comparing case A1 with cases A4 and A5, case A1 shows a slightly lower Chu's energy for both modes (10, 0) and (3, 1) upstream of x = 0.2 m. With equal initial forcings of the primary oblique and planar waves in case A1, more energy is probably transferred to accelerate the growth of secondary modes, such as the detuned one (2, 2) and subharmonic one (5, 1).

The modal growth results in figure 11 may not show straightforwardly that mode (2, 2) or (2, 0) is important in the late nonlinear stage. To provide more direct evidence, the maximum absolute contribution from mode (m, n) to the instantaneous skin friction is defined by

(4.5)\begin{equation}\Delta {C_{f,(m,n)}}(x) = \mathop {\max }\limits_{(z,t)} |{{C_{f\textrm{,(}m,n\textrm{)},disturbed}} - {C_{f,laminar}}} |,\end{equation}

where C_{f ,(m,n)},_disturbed and C_{f ,laminar} are the instantaneous skin friction induced by the laminar flow + mode (m, ±n) alone and the laminar flow alone, respectively. This indicator is defined to highlight and visualise the Fourier modes which are the source contributors to C_f in figure 9. Figure 12(a) demonstrates that mode (2, 2) is dominant in the vicinity of x = 0.3 m for both cases A1 and A2, excluding the MFD and mode (0, 2) which do not directly affect the streamwise scales. Since the large-scale $\varLambda$-vortices are observed in the same region in figure 9, it is clarified that mode (2, 2) is responsible for the formation of near-wall $\varLambda$-like structures for case A1. In addition to the baseline forcing settings, the effects of the amplitude ratio and the absolute amplitude of the initial forcings are of further interest. Figure 12(b–d) shows the results when the initial forcing ratio of the planar wave to the oblique one is set to 2:1 (case A3), 10:1 (case A4) and 1:10 (case A5), respectively. Figure 12(e) displays the result of case A6 when the absolute initial forcing decreases to 0.2 times the baseline value for both planar and oblique waves. It is observed that the pronounced region of mode (2, 2) is postponed compared with the baseline case, and that the corresponding proportion of the contribution to the skin friction becomes smaller. However, the presence of this detuned mode is found to be independent of the effects of the amplitude ratio and absolute amplitude under the considered parametric set-up.

Figure 12. Maximum absolute modal contribution to the instantaneous skin friction coefficient ΔC_f for (a) cases A1 (solid line) and A2 (dashed line), (b) case A3, (c) case A4, (d) case A5 and (e) case A6.

Figure 13 further shows a snapshot of the skin friction coefficient for cases A3–A6. Note that in cases A4 and A5, the initial planar and oblique waves are dominant, respectively. Consequently, early fundamental and oblique breakdown scenarios characterised by the aligned and staggered vortices are shown in figures 13(b) and 13(c), respectively. In terms of the combination resonance, the signature of the large-scale $\varLambda$-vortices becomes weaker if the amplitude ratio or the absolute amplitude changes. As shown in figure 13(a–d), the spanwise length scale of the $\varLambda$-vortices may be visibly affected because of other mixed-up modes, and accordingly the head angle of the vortex is altered. However, the signature of the large streamwise length scale corresponding to the detuned frequency 2ω ₀ is still maintained in the late stage of the transition, although it is postponed depending on the active region of the detuned mode (2, 2) in figure 12. These findings possibly indicate the general existence of the combination resonance mechanism if multiple instability waves are seeded. The detailed strength of this resonant scenario depends on the flow conditions, the forcing components, amplitudes, etc. The features of strength dependence and existence independence of the combination resonance have also been reported by Fezer & Kloker (2000) in supersonic and hypersonic states. In their studies, the included inflow modes (1, ±4) and (1/2, ±3) with different (relatively) low amplitudes at Mach numbers 2 and 6.8 could generate a rapidly growing ‘combination mode’ (1/2, ±7). Differently, the combination resonance under the condition of Fezer & Kloker seemed not to be pronounced enough to exceed the subharmonic and oblique breakdown scenarios somewhere and present evident flow phenomena. In addition, this mode (1/2, ±7) did not introduce new streamwise scales compared with the inflow modes. The emergence of mode (1/2, ±7) is also less complicated than that of mode (2, 2) in this paper (see (4.6)).

Figure 13. Instantaneous skin friction coefficient C_f of (a) case A3, (b) case A4, (c) case A5 and (d) case A6. The symbols λ_z,o and λ_x,o denote the spanwise and streamwise wavelengths of the oblique wave (3, 1), respectively, and λ_x,p and λ_x,(2,2) represent the streamwise wavelengths of the planar wave (10, 0) and mode (2, 2), respectively. The contour levels are identical to those in figure 9.

4.3. Weakly nonlinear stability analysis

In this section, the NPSE serves to reveal part of the physical mechanisms to explain the flow phenomena in § 4.2. As a simplified physical model for the weakly nonlinear stability analysis, only sum and difference modes of the seeded inflow disturbances are additionally generated in the NPSE simulation. Figure 14 compares the evolution of Chu's energy and profiles of the density fluctuation between the DNS and NPSE before the NPSE fails to converge. Generally, the agreement is fine. Meanwhile, the approach of the density peak between mode (3, 1) and (10, 0) supports that the evident density fluctuation of the ‘rope-like’ structure in figure 6(f) is jointly contributed to by the oblique and planar waves.

Figure 14. Comparison between DNS (solid line) and NPSE (symbol) of (a) Chu's energy and (b) profiles of density fluctuation |ρ′| at x = 0.12 m.

In addition to the self-interactions, mutual sum or difference interactions between the planar and oblique waves enable the occurrence of numerous modes. Figure 15 shows the weakly nonlinear effect on either the planar or oblique wave, where the sum or difference interaction is artificially cut off in the NPSE code. In figure 15, case B0 (see table 2) is the baseline case which is compared with DNS data. Case B0 considers the mutual interactions and self-interactions between modes (10, 0) and (3, ±1) (as well as MFD, not mentioned below). Compared with the baseline B0, case B1 interrupts the difference interaction between planar and oblique waves, which generates (7, ±1), while case B2 removes their sum interaction, which yields (13, ±1). Obviously, by comparing cases B1 and B2 with the baseline case, the difference interaction between the oblique and planar waves destabilises both oblique and planar waves, while the sum interaction stabilises both waves. The slight shift in E_Chu indicates that the direct interaction between the oblique and planar waves is weak in the initial stage. However, the two generate a series of Fourier modes via multiple generations of nonlinear interactions, which are important in the eventual breakdown. Among them, the possible detuned mode (2, 2) is of particular interest.

Figure 15. Effects of the absence of sum or difference interactions for planar waves (thick lines) and oblique waves (thin lines) on E_Chu.

Table 2. Settings for NPSE cases on the effects of sum or difference interactions.

The selection mechanism of the frequency and spanwise wavenumber for mode (2, 2) can be straightforwardly understood from the procedure of NPSE simulations. Corresponding to the dominant wavenumber 2β ₀, no other preferential frequency in the neighbourhood except for 2ω ₀ is found to be largely amplified. The detailed generation process of mode (2, 2) in the NPSE is described by

(4.6a)\begin{gather}(3,1) + (3, - 1) \to (6,0),\end{gather}

(4.6b)\begin{gather}(10,0) - (3, - 1) \to (7,1),\end{gather}

(4.6c)\begin{gather}(7,1) - (6,0) \to (1,1),\end{gather}

(4.6d)\begin{gather}(1,1) + (1,1) \to (2,2).\end{gather}

Aimed at revealing the physical mechanism of mode (2, 2), we design several NPSE cases with information given in table 3. The modes are seeded at x = 0.045 m in the NPSE with DNS-initialised profiles. Note that the primary waves have evidently larger amplitudes than the secondary instability modes in the early stage. With these primary instability waves added into the background flow, secondary instability is likely to occur. To simplify the physical problem, only primary waves (10, 0) and (3, ±1) and the concerned mode (2, 2) are included in the inflow profiles. The purpose is to identify the necessary primary instability waves which support the rapid growth of mode (2, 2).

Table 3. Settings for NPSE cases on the generation mechanism of mode (2, 2).

In table 3, cases C0 and C1 simulate the development of the small-amplitude mode (2, 2) under the large-amplitude planar and oblique waves and the base flow. Different from case C0, case C1 interrupts the sum and difference interactions between (10, 0) and (3, ±1) downstream of the inflow station x = 0.045 m. Accordingly, case C1 resembles a Floquet analysis on the secondary instability (Ng & Erlebacher 1992), where the original base flow, mode (10, 0) and mode (3, ±1) constitute a new periodic background flow. In the frame of reference moving with the group velocity of the primary waves, the combination of the primary mode (10, 0) and mode (3, ±1) varies slowly as a background flow, while small-amplitude higher-order harmonics can undergo secondary super-exponential growths. In comparison, case C2 or C3 only contains mode (10, 0) or (3, ±1) as the background primary wave, respectively. To approach the Floquet theory, a phase-locked assumption can be made. Here, we follow the usage of the terminology ‘phase-locked’ by Chang et al. (1993) in NPSE studies. Actually, ‘phase-locked’ usually refers to a state whereby two types of waves have nearly the same phase velocity (Chen et al. 2017), i.e. a lock of the phase velocity. Once the phase-locked state occurs, the generated sum or difference mode can constitute a resonant triad with the original two waves, which can give rise to a super-exponential growth. Specifically, the iterative scheme of the streamwise wavenumber except for mode (10, 0) is changed from (2.17) to

(4.7a)\begin{gather}{\mathop{\rm Im}\nolimits} ({\alpha _{mn,new}}) \!=\! {\mathop{\rm Im}\nolimits} ({\alpha _{mn,old}}) - \textrm{iRe}\left[ {\frac{1}{{{E_{mn}}}}\int_0^\infty {\bar{\rho }\left( {\hat{u}_{mn}^ \ast \frac{{\partial {{\hat{u}}_{mn}}}}{{\partial x}} \!+\! \hat{v}_{mn}^ \ast \frac{{\partial {{\hat{v}}_{mn}}}}{{\partial x}} + \hat{w}_{mn}^ \ast \frac{{\partial {{\hat{w}}_{mn}}}}{{\partial x}}} \right)\textrm{d}y} } \right],\end{gather}

(4.7b)\begin{gather}\textrm{Re}({\alpha _{mn,new}}) = \textrm{Re}({\alpha _{(10,\textrm{0})}}) \times m/10,\end{gather}

while mode (10, 0) maintains the calculation manner in (2.17). Provided that wave–wave resonance occurs, the synchronisation of the phase velocity allows the real streamwise wavenumber to be proportional to the angular frequency, as described by (4.7b).

In this complicated problem, multiple Fourier modes with approaching strength are able to affect the shape of mode (2, 2) at different wall-normal heights. Therefore, the growth of Chu's energy rather than the mode shape is chosen as the indicator to identify the dominant effect in the underlying resonant state. Figure 16 illustrates that both cases C0 and C1 can achieve an exponential growth for mode (2, 2) with a slope close to that of the DNS counterpart. The initial difference between the NPSE and DNS may come from the simplification manipulation of the NPSE in the governing equations and interaction models. Nevertheless, the agreement in the slope in the exponential amplification stage confirms that the secondary instability mechanism is mainly responsible for the growth of mode (2, 2). Furthermore, the difference between cases C0 and C1 suggests that the interaction between modes (10, 0) and (3, ±1) has a lesser effect on the amplification of mode (2, 2). The interaction between modes (10, 0) and (3, ±1) generates additional modes and constantly feeds energy to mode (2, 2) through the indirect route via these new modes in addition to a direct interaction. By comparing cases C0 or C1 with C2 and C3, it is clear that the planar or oblique primary wave alone cannot induce the secondary instability of mode (2, 2) effectively. Hence, the combination resonance of the detuned mode (2, 2) requires participation of both mode (10, 0) and mode (3, ±1). This resonant mechanism is different from conventional mechanisms, where only one primary wave is needed. We have also ruled out the rapid growth of mode (2, 2) caused only by interactions between mode (3, ±1) and other 3-D modes by additional unshown cases. In these cases, the seeded inflow modes include (3, ±1), (2, ±2) and one of the following modes to constitute an interaction triad: (1, ±1), (1, ±3), (5, ±1) or (5, ±3). The growth of mode (2, 2) in these unshown cases resembles that in cases C2 and C3. Thus, the inclusion of the primary mode (10, 0) is required to support the significant amplification of mode (2, 2). Figure 16 also shows that the phase-locked assumption essentially does not alter the results, which agrees with the conclusion by Chang et al. (1993).

Figure 16. Comparison between the NPSE results of E_Chu for cases C0 to C3 and the DNS results. Dashed line: results without the phase-locked assumption.

4.4. Resonant state and energy transfer mechanisms

The resonant state can be examined by plots of wave vectors of the Fourier modes (Fezer & Kloker 2000). If the three wave vectors (α_r, β) of the considered triad constitute a closed form, the resonance condition will be satisfied (Craik 1971):

(4.8a–c)\begin{equation}{\alpha _{r,1}} + {\alpha _{r,2}} = {\alpha _{r,3}},\quad {\beta _1} + {\beta _2} = {\beta _3},\quad {\omega _1} + {\omega _2} = {\omega _3},\end{equation}

where the angular frequency relation is automatically satisfied for the examined Fourier triad. The streamwise wavenumber is obtained via α_r = d$\varTheta$/${\rm d}\kern 0.06em x$, where $\varTheta$ is the phase angle of Fourier transform of the wall pressure in DNS. Figures 17(a) and 17(b) depict the wave vectors of the four interactions in (4.6) at x = 0.2 m and x = 0.3 m, respectively. The resonance condition for the generation of mode (2, 2) is generally fulfilled at x = 0.2 m, although the self-interaction (1, 1) + (1, 1) → (2, 2) does not fully form a closed triad. Downstream at x = 0.3 m, better closed triads indicate an enhanced resonant state. The resonance of the planar wave, the oblique wave and the detuned mode is also shown by the well-closed vector chain (2, 2) + (3, 1) + (3, −1) + (2, −2) → (10, 0) in figure 17(c).

Figure 17. Wave vectors (α_r, β) of the modes in (4.6a–d) at (a) x = 0.2 m and (b) x = 0.3 m, and (c) the vector chain (2, 2) + (3, 1) + (3, −1) + (2, −2) → (10, 0) at x = 0.3 m for case A2.

To understand the energy transfer in the nonlinear stage, a posterior energy budget analysis is presented first. With regard to mode (2, 2), the energy budget of the NPSE cases C0 to C3 in § 4.3 is considered here. Similar to the derivation by Chen et al. (2017), the wall-normal integral of the fluctuation kinetic energy equation for mode (m, n) can be written as

(4.9)\begin{equation}\int_0^\infty {\frac{1}{2}\bar{\rho }\frac{{D({{|{{u_{mn}}} |}^2} + {{|{{v_{mn}}} |}^2} + {{|{{w_{mn}}} |}^2})}}{{Dt}}\,\textrm{d}y} = \mathcal{P} + \mathcal{Q} + \mathcal{F} + \mathcal{D} + \mathcal{N},\end{equation}

where the right-hand side contains the terms in the form of a wall-normal integral, including the linear mean-shear production term $\mathcal{P}$, the non-parallel term $\mathcal{Q}$, the pressure diffusion and dilatation term $\mathcal{F}$, the viscous dissipation term $\mathcal{D}$ and the nonlinear term $\mathcal{N}$. Detailed expressions can be found in Chen et al. (2017). In addition, $\mathcal{L}$ represents the sum of all linear terms on the right-hand side.

Figure 18 indicates that the rapid energy growth of mode (2, 2) in cases C0 and C1 is caused by the positive nonlinear term. In contrast, the sums of the linear terms as well as most linear terms alone are negative and stabilise mode (2, 2), although the mean shear production effect is slight destabilisation. In terms of cases C2 and C3, the term magnitude is very small and mode (2, 2) fails to become strong.

Figure 18. Energy budget of mode (2, 2) of (a) cases C0 (solid line) and C1 (dashed line) and (b) cases C2 (solid line) and C3 (dashed line). All the terms of the four cases are normalised by the maximum $\mathcal{N}$ of case C0.

In the present study, it is difficult for the energy budget analysis to provide a comprehensive interpretation of energy transfer between multiple discrete Fourier modes. We further apply the ‘amplitude correlation method’ first proposed by Crossley et al. (1992) and Duncan & Rusbridge (1993) in plasma dynamics. Xia & Shats (2003) used this technique to provide the first experimental evidence of an inverse energy cascade in turbulence with broadband spectra. Zhang & Shi (2022) employed this method to analyse the modal interaction in the hypersonic boundary layer transition. The realisation is based on the cross-correlation function (CCF) in the form of the fourth-order covariance as

(4.10)\begin{equation}K(\tau ) = \langle [x_I^2(t)][x_{II}^2(t + \tau )]/({\parallel} [x_I^2(t)]\!\parallel{\!\parallel} [x_{II}^2(t + \tau )]\!\parallel )\rangle .\end{equation}

Here, x_I and x_II are filtered signals of the same original sample, which are centralised on the frequencies $f_{1}$ and f₂, respectively. The corresponding filters F1 and F2 are chosen as the Butterworth bandpass filter with the flattest frequency response and no phase shift. Moreover, [${\cdot}$] denotes a low-frequency filtration by filter F3 first and an operation to subtract the expected value then to obtain the fluctuation. The symbol ||${\cdot}$|| denotes the Euclidean norm for normalisation, and $\langle \cdot\rangle $ means the ensemble average. Finally, the CCF K(τ) is a function of the time lag τ. If the $\tau_{lag}$ that maximises the CCF is positive, the signal x_II lags x_I and energy flows from component f₁ to f₂ through the sideband nonlinear interaction. Otherwise, energy is considered to flow from f₂ to f₁ instead.

The physical explanation of this technique is briefly introduced here. Assume that the narrowband signals x_I and x_II contain major energy in the frequency ranges [f₁ − Δf₁, f₁ + Δf₁] and [f₂ − Δf₂, f₂ + Δf₂], respectively, similar to those in figure 11(e). Thus $x_I^2$ has two bands [2f₁ − 2Δf₁, 2f₁ + 2Δf₁] and [0, 2Δf₁] (similarly for $x_{II}^2$). The low-frequency filter F3 removes the higher band and maintains the lower band. Thus the CCF characterises information transfer between frequencies f₁ and f₂ through the sidebands with lengths of 2Δf₁ and 2Δf₂, and thus contains information on the direction of energy transfer. In the present study, the width of filter F3 is varied from 0.1f₀ to f₀ to ensure the independence of Sgn($\tau_{lag}$), where Sgn(${\cdot}$) is the sign function.

Figure 19 shows the signal samples of the wall pressure fluctuation $x_I^2$ (f₁ = 2f₀) and $x_{II}^2$ (f₂ = 19f₀) and the filtered ones by F3 as well as the CCF. The sampled station is at x = 300 mm, where the strong nonlinear interaction gives rise to stochastic time signals and the concerned $\varLambda$-vortices are evident. The pressure signal is chosen for the analysis to approach an experimental signal (Zhang & Shi 2022). The two example results of the CCFs in figure 19(c) indicate that the frequency f₁ = 2f₀ lags both frequencies f₂ = 19f₀ and f₂ = 10f₀, and thus energy flows from these two high frequencies to f₁ = 2f₀.

Figure 19. Signal sample of wall pressure fluctuations at x = 300 mm for case A2 of (a) $x_I^2$ for f = 2f₀ and (b) $x_{II}^2$ for f = 19f₀, and (c) the resulting CCF. Thick lines in (a,b): filtered signal by the low-frequency filter F3. Dashed-dotted line in (c): another CCF plot with $x_{II}^2$ belonging to f = 10f₀.

By means of the approach above, how mode (2, 2) sustains its energy and supports the formation of $\varLambda$-vortices can be analysed. Let f₁ = 2f₀, and the characteristic time lag $\tau_{lag}$ is obtained via the calculation of CCF for each pair of (f₁, f₂). If Sgn($\tau_{lag}$) > 0, energy flows from f₁ to f₂; if Sgn($\tau_{lag}$) < 0, energy flows from f₂ to f₁; if Sgn($\tau_{lag}$) = 0, there is no evident peak of the CCF with a non-negligible cross-correlation, say, (CCF)_max < 0.1. Figure 20 shows a plot of Sgn($\tau_{lag}$) and the resulting diagram of the energy flow direction. It is known from the energy budget analysis in figure 18 that mean shear slightly destabilises mode (2, 2) and other linear terms stabilise it. Therefore, mean shear feeds energy to mode (2, 2), and mode (2, 2) dissipates energy into other linear mechanisms. Meanwhile, the constantly growing MFD is the product of the sum interaction between each mode and its complex conjugate, i.e. (2, 2) + (−2, −2) → (0, 0). It is thus inferred that mode (2, 2) participates in energy transfer to the MFD. The remaining groups in figure 20(b) are categorised into nonlinear energy ‘contributors’ and ‘receivers’. In the studied case, the energy contributors for f₁ = 2f₀ consist of two neighbours f₂ = f₀ and f₂ = 4f₀, the optimal planar wave frequency f₂ = 10f₀ and additional higher modes from f₂ = 12f₀ to f₂ = 20f₀. Note that the low frequency f₂ = f₀ is also one of the source modes for mode (2, 2) via (1, 1) + (1, 1) → (2, 2) (see (4.6) for details). The energy receivers for mode (2, 2) are the optimal oblique wave frequency f₂ = 3f₀ and other components f₂ = 8f₀, f₂ = 9f₀ and f₂ = 11f₀ in the vicinity of the frequency of mode (10, 0). Corresponding to each energy flow direction from f₂ to f₁ or from f₁ to f₂, there are two types of potential interaction processes to allow nonlinear energy transfer: the sum interaction (f₁ + f₂ → f₃ = f₁ + f₂) and difference interaction (f₁ − f₂ → f₃ = |f₁ − f₂|) in conjunction with their reverse counterparts, such as (f₁ + f₂) − f₁ → f₂. In summary, figure 20(b) gives a physical image of energy transfer through linear mechanisms and through the narrow sidebands in nonlinear modal interactions.

Figure 20. The energy flow direction at x = 0.3 m for f₁ = 2f₀ of case A2: (a) sign function of the time lag of f₂ compared to f₁ and (b) diagram of the energy flow direction.

Note that the energy flow direction may depend on the signal type and wall-normal height for analysis. Amplitude correlation analysis is further performed based on the (ρu)′ signal at the height of its maximum for mode (3, 1) at x = 300 mm. What is different is that f₁ = f₀ and the high frequencies f₁ = 13f₀, 14f₀ and 15f₀ change from ‘contributors’ to ‘receivers’, while f₁ = 3f₀ changes from ‘receivers’ to ‘contributors’. However, the behaviour that multiple frequencies participate in nonlinear energy transfer as ‘contributors’ and ‘receivers’ is found to be consistent.

5. Conclusions

Concentrating on the linear and nonlinear dynamics of optimal disturbances, a Mach 6 flat-plate boundary layer is studied by resolvent analysis, PSE and DNS. Newly observed flow phenomena and physical interpretations of the modal interaction and breakdown scenario are presented.

The resolvent analysis demonstrates that the optimal planar wave mode (10, 0) and the optimal oblique wave mode (3, 1) are notably more significant than the optimal streak. Due to the strong transient growth, the integrated Chu's energy of the oblique wave is comparable with that of the planar wave, although the growth rate of the second mode exceeds that of the first mode. Introducing the two optimal disturbances jointly with the same magnitude of initial forcing, the DNS results show that the oblique wave is highly pronounced in the outer boundary layer. Quantitatively, the active performance of the non-acoustic oblique wave away from the wall may not be fully detected by conventional wall pressure measurements, such as PCB sensors.

Modal decomposition via a Fourier transform identifies the coexistence of the oblique breakdown, fundamental resonance triggered by mode (10, 0), combination resonance related to a detuned mode (2, 2) and subharmonic resonance triggered by mode (10, 0). Particularly, mode (2, 2) is responsible for the occurrence of large-scale staggered $\varLambda$-vortices near the wall before the eventual breakdown. More importantly, the presence of mode (2, 2) and the $\varLambda$-vortices are found not to depend on the (considered) amplitude ratio and absolute amplitude of the initial forcings corresponding to the oblique and planar waves. This suggests that the combination resonance can be a viable transition scenario with the presence of multiple instability waves. The NPSE analyses further confirm that the detuned mode (2, 2) is mainly a consequence of the secondary instability, where the background primary waves require participation of both modes (10, 0) and (3, 1). The nonlinear interaction between the two types of primary waves generates a series of discrete Fourier modes. These multiple components of spectra in the early nonlinear stage are different from the observation in a single common breakdown scenario. In addition, the resonant state associated with these multiple components is observed. Finally, the energy budget and amplitude correction techniques describe the energy flow direction. The linear production and dissipation effects as well as nonlinear interactions with ‘contributor’ modes and ‘receiver’ modes are involved in the physical image of energy transfer for mode (2, 2). The present study indicates that the boundary layer transition induced by multiple primary instability waves can contain rich flow physics and unexpected breakdown scenarios. Particularly, the significance of combination resonance and primary oblique wave in a hypersonic state is highlighted in this work. The multiple resonance mechanisms including the combination one, revealed by a specific initialisation of frequency and wavenumber pair in the present paper, may constitute a building block in a real natural transition with a wide and continuous frequency range.