2.1. Direct numerical simulation

The governing equation for the computation of flow fields is the compressible Navier–Stokes equation in the conservation form

(2.1) \begin{align} \frac {{\partial {\boldsymbol{U}}}}{{\partial t}} + \frac {{\partial {\boldsymbol{F}}}}{{\partial x}} + \frac {{\partial {\boldsymbol{G}}}}{{\partial y}} + \frac {{\partial {\boldsymbol{H}}}}{{\partial z}} = {\unicode{x1D63D}\boldsymbol{f'}}, \end{align}

where ${\boldsymbol{U}} = (\rho , \,\rho u,\,\rho v,\,\rho w,\,\rho e)^{\textrm {T}}$ is the vector of conservative variables, and $\boldsymbol{F}$ , $\boldsymbol{G}$ and $\boldsymbol{H}$ are the vectors of inviscid and viscous fluxes. Detailed expressions can be found in Anderson (Reference Anderson1995). Here, $\rho$ is density, and $u$ , $v$ and $w$ are velocities in the Cartesian $x$ , $y$ and $z$ directions, respectively. Total energy per unit mass is denoted by $e$ . The matrix $\unicode{x1D63D}$ constrains the forcing vector ( $\boldsymbol{f'}$ ) to be added at a certain region, which is set to be $x$ = 0.04 m in this study. In the DNS, the optimal forcings obtained by resolvent analysis are employed to initiate the unsteady flow.

The perfect gas model is employed with a constant specific heat ratio of 1.4. Furthermore, the dynamic viscosity is calculated using Sutherland’s law with a constant $T_s$ = 110.4 K, and the thermal conductivity coefficient is calculated based on a constant Prandtl number 0.72. The simulation constitutes two steps. In the first step, the base flow is calculated, and the right-hand side term of (2.1) is set to zero. Subsequently, the resolvent analysis or DNS is performed on the converged base flow. It is assumed that the vector of conservation variable can be decomposed into its base-flow and perturbation parts

(2.2) \begin{align} {\boldsymbol{U}}(x,y,z,t) = {\boldsymbol{\bar U}}(x,y) + {\boldsymbol{U'}}(x,y,z,t). \end{align}

Hereafter, the overbar and prime represent mean variables and perturbation variables, respectively.

The DNS is performed using a multi-block parallel finite difference solver called OpenCFD, which has been successfully employed in high-speed transitional and turbulent simulations (Li et al. 2002, Li et al. 2009). The inviscid flux splitting is implemented by the Steger – Warming scheme, and a seven-order weight essential non-oscillation scheme is employed for the reconstruction of the split flux. The sixth-order central difference scheme is utilised for viscous flux discretisation, and a third-order Runge – Kutta method is employed for temporal integration. With regard to the boundary condition, the left boundary of the computational domain is free stream, and the wall is set to be isothermal, no slip and no penetration or with TDIBC in a certain region. Furthermore, the computational domain extends upstream to $x$ = −0.05 m with a slip-wall (inviscid symmetry) boundary condition to ensure proper inflow development and mitigate potential numerical instabilities. The upper and outer boundaries are addressed by extrapolation. The spanwise boundary is set to be periodic.

2.2. Time-domain impedance boundary condition

Direct numerical simulation with realistic metasurface microstructures can be conducted to investigate the transition control mechanism. However, it is computationally expensive to mesh numerous micro-cavities. To efficiently and simultaneously simulate the interactions between acoustic metasurfaces and multi-modal waves, the application of an IBC is recommended. The early acoustic impedance models of a metasurface were all built in the frequency domain and applied with a single disturbance frequency in hypersonic boundary-layer transition studies (Fedorov et al. Reference Fedorov, Malmuth and Rasheed2001; Möser & Michael Reference Möser2009), which means that they are not suitable for transition considering broadband disturbances. In this study, a multi-pole broadband impedance model with a piecewise linear recursive convolution technique proposed by Fung & Ju (Reference Fung and Ju2004) will be revised and incorporated into the current DNS code afterwards. This embedded boundary condition enables efficient simulations of broadband-disturbance propagation by transforming the acoustic IBC from the frequency domain into the time domain.

The reflective coefficient $\tilde R$ in the frequency domain is defined as

(2.3) \begin{align} \tilde R(\omega ) = {{{A^{\textit {out}}}(\omega )} \mathord {\left / {\vphantom {{{A^{\textit {out}}}(\omega )} {{A^{\textit {in}}}}}} \right . } {{A^{\textit {in}}}}}(\omega ), \end{align}

where $\omega$ is the angular frequency, ${A^{\textit {in}}}$ and ${A^{\textit {out}}}$ are the amplitudes of ingoing and outgoing waves, respectively. The ingoing and outgoing waves are related to the physical variables of the flow field (pressure and velocity) at the wall (Fung & Ju Reference Fung and Ju2004). Their amplitudes are defined as

(2.4) \begin{align} {A^{\textit {in}}}(t) = v'(t) - p'(t)\textit{Ma}_{\infty }\sqrt {{T_{w}}} ,\,{A^{\textit {out}}}(t)= v'(t) + p'(t)\textit{Ma}_{\infty }\sqrt {{T_{w}}}. \end{align}

where $\textit{Ma}$ is the Mach number and $p$ is the pressure. The wall softness $\tilde S$ is given by

(2.5) \begin{align} \tilde S(\omega ) = 1 + \tilde R(\omega ). \end{align}

Physically $\tilde S$ = 2 and $\tilde S$ = 1 when the wall is totally reflective (rigid wall) and absorptive (soft wall), respectively. Combining (2.3) and (2.5) gives rise to

(2.6) \begin{align} {A^{\textit {out}}}(\omega ) = - {A^{\textit {in}}}(\omega ) + \tilde S(\omega ){A^{\textit {in}}}(\omega ). \end{align}

The softness can be approximated by a multi-oscillator model proposed by Fung & Ju (Reference Fung and Ju2004)

(2.7) \begin{align} \tilde S(s) = \mathop \sum \limits _{k = 1}^{{n_0}} \left [ {\frac {{{\mu _k}}}{{s - {p_k}}} + \frac {{\mu _k^\dagger }}{{s - p_k^\dagger }}} \right ] + {C_0},\: s = {\textrm {i}}\omega , \end{align}

where $p_k$ and $p_k^\dagger$ are poles of pole base function, and ${\mu _k} = {\textrm {i}} \boldsymbol{\cdot }{\textrm {Residue}} [ {\tilde S(s),{p_k}} ]$ . The superscript $\dagger$ refers to complex conjugate, and $n_0$ is the total number of pole pairs. In addition, $C_0$ is a constant real number. The utilisation of this constant number to improve the accuracy of the model was also reported by Fiévet et al. (Reference Fiévet, Deniau, Brazier and Piot2021). Due to the requirement of the reality of the signal (Rienstra Reference Rienstra2006), the softness in the time domain must be real and ${\tilde S^\dagger }(s) = \tilde S( - s)$ . The values of poles and residue are obtained from a nonlinear fitting process (Douasbin et al. Reference Douasbin, Scalo, Selle and Poinsot2018) to approximate the softness calculated from a well-validated impedance model, which considers the effect of high-order diffraction (Zhao et al. Reference Zhao, Liu, Wen and Cheng2018a ). The authors have fully validated the impedance model by comparing it with the DNS result over meshed metasurface (Zhao, Wen & Long Reference Zhao, Wen, Long, Tian, Zhou and Wu2019; Guo et al. Reference Guo, Liu, Zhao, Hao and Wen2023b ). Then, the softness in the time domain can be obtained using the inverse Fourier transform

(2.8) \begin{align} \tilde S(t) = \frac {1}{{2\pi }}\int _{ - \infty }^\infty {\tilde S(\omega ){e^{{\textrm {i}}\omega t}}} {\rm d}\omega . \end{align}

Using the Residue theorem, one may obtain

(2.9) \begin{align} \tilde S(t) = \left ( {\mathop \sum \limits _{k = 1}^{{n_0}} {\mu _k}{e^{{p_k}t}} + \mathop \sum \limits _{k = 1}^{{n_0}} \mu _k^\dagger {e^{p_k^\dagger t}}} \right )H(t) + {C_0}\delta (t). \end{align}

Here, the Heaviside function $H$ ( $t$ ) indicates that the causality condition is satisfied, and $\delta$ is the Dirac delta function. Taking the inverse Fourier transform of (2.6), one may obtain

(2.10) \begin{align} {A^{\textit {out}}}(t) = ({C_0} - 1){A^{in}}(t) + \int _0^\infty {\tilde S(\tau ){A^{in}}(t - \tau )}{\rm d} \tau . \end{align}

where $\tau$ is the dummy integration variable representing time. Substituting (2.9) into (2.10) yields

(2.11) \begin{align} {A^{\textit {out}}}(t) = ({C_0} - 1){A^{\textit {in}}}(t) + \int _0^\infty {\mathop \sum \limits _{k = 1}^{{n_0}} \left({\mu _k}{e^{{p_k}\tau }} + \mu _k^\dagger {e^{p_k^\dagger \tau }} \right){A^{\textit {in}}}(t - \tau )} {\rm d}\tau . \end{align}

For the next time at step $t + \Delta t$ , (2.11) can be written as

(2.12) \begin{align} {A^{\textit {out}}}(t + \Delta t) = ({C_0} - 1){A^{\textit {in}}}(t + \Delta t) + \int _0^\infty {\mathop \sum \limits _{k = 1}^{{n_0}} \left({\mu _k}{e^{{p_k}\tau }} + \mu _k^\dagger {e^{p_k^\dagger \tau }}\right){A^{\textit {in}}}(t + \Delta t - \tau )} {\rm d}\tau . \end{align}

Here, the convolution in (2.12) can be written as $G_k^{\textit {in}}$ , that is

(2.13) \begin{align} G_k^{\textit {in}}(t + \Delta t) = \int _0^\infty {{\mu _k}{e^{{p_k}\tau }}{A^{\textit {in}}}(t + \Delta t - \tau )} {\rm d}\tau . \end{align}

According to Fung & Ju (Reference Fung and Ju2004), (2.13) can be calculated by a recursive scheme

(2.14) \begin{align} G_k^{\textit {in}}(t + \Delta t) = {z_k}G_k^{\textit {in}}(t) + {\mu _k}\Delta t \big[{w_{k0}}{A^{\textit {in}}}(t + \Delta t) + {w_{k1}}{A^{\textit {in}}}(t) \big]. \end{align}

Here, ${z_k} = {e^{{p_k}\Delta t}}$ , and

(2.15) \begin{align} \left \{ \begin{aligned} w_{k0} &= \frac {z_k - 1}{p_k^2 \Delta t^2} - \frac {1}{p_k \Delta t}, \\ w_{k1} &= -\frac {z_k - 1}{p_k^2 \Delta t^2} + \frac {z_k}{p_k \Delta t}. \end{aligned} \right . \end{align}

With a first-order approximation (Scalo et al. Reference Scalo, Bodart and Lele2015), the ingoing wave at $t + \Delta t$ can be expressed by

(2.16) \begin{align} {A^{\textit {in}}}(t + \Delta t) \simeq {A^{\textit {in}}}(t) + \sqrt {{T_w}} / \textit{Ma}_{\infty }\Delta t\frac {\partial }{{\partial y}}\big ( { v'(t) - \textit{Ma}_{\infty }\sqrt {{T_w}} p'(t)} \big ). \end{align}

Then, combining (2.12) and (2.16) gives rise to the fluctuating wall-normal velocity at $t + \Delta t$

(2.17) \begin{align} v'(t + \Delta t) & = \big ( {{A^{\textit {in}}}(t + \Delta t) + {A^{out}}(t + \Delta t)} \big )/2 \nonumber \\ & = {\textrm {Re}}\left ( {\mathop \sum \limits _{k = 1}^{{n_0}} G_k^{\textit {in}}(t + \Delta t)} \right ) + {C_0}{A^{\textit {in}}}(t + \Delta t). \end{align}

Finally, the pressure fluctuation at $t + \Delta t$ is

(2.18) \begin{align} p'(t + \Delta t) = {{\big [ { v'(t + \Delta t)-{A^{\textit {in}}}(t + \Delta t)} \big ]} \mathord {\left / {\vphantom {{\left [ {{A^{in}}(t + \Delta t) + v'(t + \Delta t)} \right ]} {\left ( {\textit{Ma}_{\infty }\sqrt {{T_{T_{w}}}} } \right )}}} \right . } {\big ( {\textit{Ma}_{\infty }\sqrt {{T_w}} } \big )}}. \end{align}

In the simulation where the metasurface exists locally, the pressure and normal velocity on the wall are updated based on (2.17) and (2.18) at each time step, and the incoming wave amplitude ${A^{\textit {in}}}(t)$ is initialised using (2.4). The density is calculated by the updated pressure from the equation of state for perfect gases. Although the TDIBC was derived from linear approximations, it remains valid for the nonlinear disturbance dynamics in our system due to the strongly viscous-dominated regime inside the porous medium (Chen & Scalo Reference Chen and Scalo2021a ; Sousa et al. Reference Sousa, Wartemann and Wagner2024).

2.3. Resolvent analysis

The resolvent analysis provides the most energetic response of the flow field due to per unit energy of external forcing. Substituting (2.2) into (2.1), and subtracting the base-flow equation yield the following form:

(2.19) \begin{align} \frac {{\partial {\boldsymbol{U'}}}}{{\partial t}} + \frac {{\partial {\boldsymbol{F'}}}}{{\partial x}} + \frac {{\partial {\boldsymbol{G'}}}}{{\partial y}} + \frac {{\partial {\boldsymbol{H'}}}}{{\partial z}} = {\unicode{x1D63D}\boldsymbol{f'}}. \end{align}

Considering a small-amplitude forcing term to study the instability, one may obtain

(2.20) \begin{align} \frac {{\partial {\boldsymbol{U'}}}}{{\partial t}} = {\boldsymbol{{\unicode{x1D63C}}U'}} + {\unicode{x1D63D}\boldsymbol{f'}}, \end{align}

where matrix $\unicode{x1D63C}$ is the Jacobian matrix constituted by the base-flow variables. The harmonic assumption is utilised for the small-amplitude perturbation vector

(2.21) \begin{align} {\boldsymbol{U'}}(x,y,z,t) = {\boldsymbol{\hat U}}(x,y)\exp ({\textrm {i}}\beta z - {\textrm {i}}\omega t) + {\textrm {c}}{\textrm {.c}}., \end{align}

where $\boldsymbol{\hat U}$ is the complex eigenfunction, $\beta$ is the spanwise wavenumber, $\omega$ is the angular frequency and ${\textrm {c}}{\textrm {.c}}.$ denotes complex conjugate. Similarly, the harmonic forcing can be written as

(2.22) \begin{align} {\boldsymbol{f'}}(x,y,z,t) = {\boldsymbol{\hat f}}(x,y)\exp ({\textrm {i}}\beta z - {\textrm {i}}\omega t) + {\textrm {c}}{\textrm {.c}}., \end{align}

Substituting (2.21) and (2.22) into (2.20) gives

(2.23) \begin{align} {\boldsymbol{\hat U}} = {\unicode{x1D64D}}{\unicode{x1D63D}}\boldsymbol{\hat f},{\begin{array}{*{20}{c}} {}&{{\unicode{x1D64D}} = ( - {\textrm {i}}\omega {\unicode{x1D644}} - {\unicode{x1D63C}})} \end{array}^{\!\!\!-1}}, \end{align}

where $\unicode{x1D64D}$ represents the response matrix, indicating the relationship between the external forcing and the linear response of the system. Here, the identity operator is represented by $\unicode{x1D644}$ . In resolvent analysis, the maximal amplification of the energy, i.e. the optimal gain $\sigma ^2$ , is sought in the parameter space of ( $\omega$ , $\beta$ ). The optimal gain is defined by the energy ratio of the output (response) to the input (forcing) of the system

(2.24) \begin{align} {\sigma ^2}(\beta ,\omega ) = \mathop {\max }\limits _{{\boldsymbol{\hat f}}} \frac {{{{\big \| {{\boldsymbol{\hat U}}} \big \|}_E}}}{{{{\big \| {{\unicode{x1D63D}\boldsymbol{\hat f}}} \big \|}_E}}}. \end{align}

In this study, Chu’s energy (Chu & Kovasznay Reference Chu and Kovásznay1958) is employed to calculate the energy norm

(2.25) \begin{align} {\big \| {{\boldsymbol{\hat U}}} \big \|_E} = \iint \limits _{\varOmega } {{{\boldsymbol{U}}^\bot }{\unicode{x1D648}\boldsymbol{\hat U}}}{\rm d}x{\rm d}y, \end{align}

where $\varOmega$ represents the computational domain for resolvent analysis, the superscript $ \bot$ refers to the conjugate transpose and $\unicode{x1D648}$ is the weight operator given by Bugeat et al. (Reference Bugeat, Chassaing, Robinet and Sagaut2019). A local energy integral in a dimensionless form is defined by

(2.26) \begin{align} {E_{\textit{Chu}}}(x) = \frac {1}{2}\int _0^\infty {\left [ {\bar \rho \big({{u'}^2} + {{v'}^2} + {{w'}^2} \big) + \frac {{\bar T}}{{\gamma \textit{Ma}_\infty ^2\bar \rho }}{{\rho '}^2} + \frac {{\bar \rho }}{{\gamma (\gamma - 1)Ma_\infty ^2\bar T}}{{T'}^2}} \right ]} \,{\rm d}y. \end{align}

The optimisation problem in (2.24) can be transformed into an eigenvalue problem with respect to $\sigma ^2$ , as demonstrated by Sipp & Marquet (Reference Sipp and Marquet2013). The resulting discrete eigenvalue problem is solved using ARPACK software for given values of $\beta$ and $\omega$ associated with the regular mode (Sorensen et al. Reference Sorensen, Lehoucq, Yang and Maschhoff1996). Additional details regarding the resolvent analysis solver and the associated validation cases can be found in Hao et al. (Reference Hao, Cao and Guo2023) and Guo et al. (Reference Guo, Liu, Zhao, Hao and Wen2023b ).

2.4. Linear stability theory and parabolised stability equation

To identify transient growth and modal growth captured by resolvent analysis, LST and parabolised stability equation (PSE) are employed. The PSE further considers the non-parallel effect compared with LST. Specifically, the LST provides the initial inlet profiles (eigenfunctions) for the PSE, and the PSE can obtain the non-local evolution of eigenmodes, including the first and second modes of interest here. In PSE, the disturbance $\boldsymbol{\psi '}$ is expressed by

(2.27) \begin{align} {\boldsymbol{\psi '}}(x,y,z,t) = {\boldsymbol{\hat \psi }}(x,y)\exp \left ( {{\textrm {i}}\int _{{x_0}}^x {\alpha {\rm d}\tilde x} + {\textrm {i}}\beta z - {\textrm {i}}\omega t} \right )+ {\textrm {c}}{\textrm {.c}}., \end{align}

where the vector ${\boldsymbol{\psi }} = {(\rho ,u,v,w,T)^{\textrm {T}}}$ , $\boldsymbol{\hat \psi }$ and $\alpha$ are the shape function and the complex streamwise wavenumber, respectively, and $x_0$ is the initialisation location of PSE marching. Substituting (2.27) into (2.19) and dropping the forcing terms give rise to the PSE governing equation

(2.28) \begin{align} ({{\unicode{x1D647}}_0} + {{\unicode{x1D647}}_1}){\boldsymbol{\hat \psi }} + {{\unicode{x1D647}}_2}\frac {{\partial {\boldsymbol{\hat \psi }}}}{{\partial x}} + \frac {{\partial \alpha }}{{\partial x}}{{\unicode{x1D647}}_3}{\boldsymbol{\hat \psi }} = \boldsymbol{0}. \end{align}

The effects of the locally parallel flow, the non-parallel base flow, the non-local shape function and the streamwise-varying wavenumber are absorbed in the base-flow-related operators ${\unicode{x1D647}}_0$ , ${\unicode{x1D647}}_1$ , ${\unicode{x1D647}}_2$ and ${\unicode{x1D647}}_3$ , respectively. Detailed expressions of these operators can be found in Paredes (Reference Paredes2014). An eigenvalue problem is solved in LST when keeping only the local operator ${\unicode{x1D647}}_0$ in (2.28). The calculation is performed by an in-house code CHASES, which integrates the LST, linear parabolised stability equation (LPSE), nonlinear parabolised stability equation and sensitivity analyses. The code has been validated by a series of cases for LST and LPSE compared with both theoretical (Guo et al. Reference Guo, Gao and Jiang2020, Guo et al. Reference Guo, Gao and Jiang2021, Guo et al. Reference Guo, Shi, Gao, Jiang, Lee and Wen2022b) and DNS (Cao, Hao & Guo Reference Cao, Hao, Guo, Wen and Klioutchnikov2023; Hao et al. Reference Hao, Cao and Guo2023; Guo et al. Reference Guo, Hao and Wen2023a ) results. The detailed formulation and the numerical method can be found in the references (Paredes Reference Paredes2014; Guo et al. Reference Guo, Hao and Wen2023a ).

2.5. Case initialisation for direct numerical simulation

To explore the effect of the acoustic metasurface in hypersonic transitional flows, three cases subject to the same forcing yet different wall boundary conditions are considered here. In the present Mach 6 state, Guo et al. (Reference Guo, Hao and Wen2023a ) reported two dominant optimal disturbances: a low-frequency oblique wave and a high-frequency planar wave. These two types of forcing are considered in the transition simulation. The mathematical form of the forcing in DNS is given by

(2.29) \begin{align} {\boldsymbol{f'}}({x_0},y,z,t) = {\boldsymbol{f'}_{\!p}}({x_0},y,z,t) + {\boldsymbol{f'}_{\!o}}({x_0},y,z,t), \end{align}

where $x_0$ = 0.04 m is the specified forcing location, and subscripts ‘ $p$ ’ and ‘ $o$ ’ represent the optimal planar wave and oblique wave obtained by resolvent analysis, respectively. According to our previous study (Guo et al. Reference Guo, Hao and Wen2023a ), the additional background noise term produces a neglectable effect on the flow field, which is thus not considered here. The planar-wave forcing is detailed as

(2.30) \begin{align} {\boldsymbol{f'}_{\!p}}({x_0},y,z,t) = {\varepsilon _{_p}}{{{\boldsymbol{\hat f}}}_{\!p}}({x_0},y)\exp ({\textrm {i}}{\beta _p}z - {\textrm {i}}{\omega _p}t) + {\textrm {c}}{\textrm {.c}}., \end{align}

where $\varepsilon$ is the amplitude coefficient, and $\beta _p = 0$ for the planar wave employed . A pair of oblique waves with an opposite spanwise wavenumber is employed, which is expressed by

(2.31) \begin{align} {\boldsymbol{f'}_{\!o}}({x_0},y,z,t) = {\varepsilon _{_o}}\big [ {{{{\boldsymbol{\hat f}}}_o}({x_0},y)\exp ({\textrm {i}}{\beta _o}z - {\textrm {i}}{\omega _o}t) + {{{\boldsymbol{\hat f}}}_{\!-o}}({x_0},y)\exp ( - {\textrm {i}}{\beta _o}z - {\textrm {i}}{\omega _o}t)} \big ] + {\textrm {c}}{\textrm {.c}}, \end{align}

In all transition simulation cases, the initial amplitude $(\rho u)'_{\textit{max}}$ = 0.002 remains the same for planar wave $({\omega _p},{\beta _p})$ and a pair of oblique waves $({\omega _o}, \pm {\beta _o})$ at $x$ = 0.045 m. This set-up excites the first mode and second mode of equal importance near the forcing. In this case, strong mutual interactions can occur with the presence of two primary instabilities. As reported by Guo et al. (Reference Guo, Hao and Wen2023a ), the resulting nonlinear mechanism (combination resonance) is independent of both the initial amplitude ratio and the absolute amplitude of the oblique and planar waves. Therefore, this study focuses on one specified set-up of wave amplitudes without loss of generality.

The wall boundary conditions for DNS cases are presented in table 1. The boundary condition is of primary interest in this study. In case 1, only the solid-wall boundary condition is employed. For case 2 with the acoustic metasurface, TDIBC is employed starting from $x$ = 0.115 m to the end of the computational domain, i.e. $x$ = 0.6 m. It is important to note that $x = 0.115$ m corresponds to the neutral point of the optimal plane wave (10, 0), downstream of which the modal growth begins. In case 3, during the late stage of the transitional flow ( $x$ >0.34m), the TDIBC is replaced by solid walls to mitigate undesirable higher wall friction and flux in the late transitional stage. The selection of $x$ = 0.34 m in case 3 also considers the modal growth region of the second mode. As illustrated in figure 2 ( $b$ ), the modal growth of the optimal plane wave (10, 0) terminates at approximately $x$ = 0.34 m. A detailed discussion and explanation will be provided in the subsequent sections.

Table 1. Streamwise range of different wall boundary conditions for transitional DNS cases.

In the DNS, the computational domain ranges from $x$ = 0 to $x$ = 0.6 m. The number of the grid points is $3060 \times 260 \times 60$ , denoting the grid numbers in the streamwise, wall-normal and spanwise directions, respectively. The corresponding dimensionless grid spacings are $\Delta {x^ + \approx }$ 3.13, $\Delta z^ + \approx$ 5.9 and $\Delta y_{\textit{min}}^ + \approx$ 0.30, which are evaluated based on the procedure described by Guo et al. (Reference Guo, Shi, Gao, Jiang, Lee and Wen2022a , Reference Guo, Hao and Wen2023a ). Grid convergence of the amplitude evolution of the optimal disturbances has been confirmed (Guo et al. Reference Guo, Hao and Wen2023a ). A good agreement is also reached between the DNS and the stability analysis.

2.6. Bispectral mode decomposition

Based on the DNS data, bi-Fourier analysis and bispectral analysis are performed. Bi-Fourier analysis demonstrates the significant amplitude evolution for multiple disturbances, while bi-spectral decomposition reveals the energy-transfer mechanism. The BMD method developed by Schmidt (Reference Schmidt2020) is employed to elucidate the triadic interactions arising from quadratic nonlinearities in the Navier–Stokes equations, which govern energy transfer between scales in the flow. As shown in Schmidt (Reference Schmidt2020), BMD examines triadic interactions originating from the quadratic nonlinearity inherent in the Navier–Stokes equations. The perturbation dynamics about a base flow is governed by

(2.32) \begin{align} \frac {\partial \boldsymbol{q}'}{\partial t} = \unicode{x1D647}\boldsymbol{q}' + \unicode{x1D64C} \big(\boldsymbol{q}',\boldsymbol{q}'\big), \end{align}

where $\boldsymbol{q}'(\boldsymbol{x},t)$ represents the flow perturbation (typically velocity components or pressure fluctuations), $\unicode{x1D647}$ is the linear operator and $\unicode{x1D64C}(\boldsymbol{\cdot },\boldsymbol{\cdot })$ denotes the quadratic nonlinearity arising from terms like $(\boldsymbol{u}' \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}'$ in the convective term.

The expectation operator $E[\boldsymbol{\cdot }]$ , central to BMD’s statistical formulation, performs ensemble averaging over independent flow realisations. Specifically, BMD identifies phase-coupled frequency triads $\{f_k, f_l, f_m\}$ satisfying $f_m = f_k \pm f_l$ through the bispectrum. The bispectrum is defined as the double Fourier transform of the third-order moment $E[q(t)q(t-{\tau _k})q(t-{\tau _l})]$

(2.33) \begin{align} S_{\textit{qqq}}(f_k, f_l) = \lim _{\mathcal{T} \to \infty } \frac {1}{\mathcal{T}} E \big [ \hat {q}^{\dagger }(f_k) \hat {q}^{\dagger }(f_l) \hat {q}(f_k + f_l) \big ], \end{align}

where $\hat {\boldsymbol{q}}(\boldsymbol{x},f)$ is the Fourier transform of ${\boldsymbol{q}}(\boldsymbol{x},f)$ and $\dagger$ indicates complex conjugate. The classical bispectrum $S_{\textit{qqq}}$ serves as a conceptual foundation for detecting quadratic phase coupling. Its spatial generalisation, the mode bispectrum $\phi$ , is achieved by replacing the temporal expectation $E[\boldsymbol{\cdot }]$ with a spatially weighted ensemble average over Fourier realisations, thereby extending pointwise statistics to coherent structures. For multidimensional datasets, BMD constructs the bispectral density tensor

(2.34) \begin{align} \unicode{x1D63F} = \frac {1}{N_{\textit{blk}}} \hat {\unicode{x1D64C}}_{k \circ l}^H \unicode{x1D652} \hat {\unicode{x1D64C}}_{k+l}, \end{align}

with $\unicode{x1D652}$ containing spatial quadrature weights, $N_{\textit{blk}}$ the number of data blocks and $\hat {\unicode{x1D64C}}_{k \circ l} = \hat {\unicode{x1D64C}}_k \circ \hat {\unicode{x1D64C}}_l$ being the Hadamard product of Fourier realisations. Here, superscript $H$ denotes the Hermitian transpose (conjugate transpose). The matrix $\unicode{x1D64C}$ organises the Fourier coefficients of $\boldsymbol{q}'(\boldsymbol{x},t)$ in a block form, where each column contains a realisation of $\hat {\boldsymbol{q}}(\boldsymbol{x},f)$ for spectral estimation. The dominant triadic interactions are identified by maximising the numerical radius of $\unicode{x1D63F}$ , which involves solving the eigenvalue problem for the Hermitian matrix $\unicode{x1D643}(\theta ) = {1}/{2}(e^{i\theta } \unicode{x1D63F} + e^{-i\theta } \unicode{x1D63F}^H)$ .The complex mode bispectrum $\phi (f_k, f_l)$ , which is defined as the maximum numerical radius of $\unicode{x1D63F}$ , quantifies the triad interaction strength through its magnitude $|\phi |$ . The magnitude of the mode bispectrum scales with the product of the absolute amplitudes of the triad components. This preserves the energy-transfer scaling inherent to quadratic nonlinearities.

The BMD method and open-source code proposed by Schmidt (Reference Schmidt2020) were employed in this study, and a detailed introduction for the formulation and definition can be found in Schmidt (Reference Schmidt2020).

2.7. A summary of analytical techniques

The present analysis framework systematically combines existing techniques to elucidate the effects of the metasurface across different transition stages. The investigation progresses from linear mechanisms, in which the resolvent analysis identifies optimal disturbances and LST/PSE quantifies modal growth rates, to the nonlinear transition dynamics captured via DNS visualisations of vortex formation and bi-Fourier quantifications of mode evolutions with or without the metasurface. The complete energy pathway is resolved through BMD of triadic interactions and energy budget diagnosis of dissipation mechanisms. Crucially, DNS visualisations provide physical grounding for spectral analyses, with instantaneous flow structures associated with BMD-detected triadic interactions and streak spacing matching bi-Fourier analysis. This represents the first coupled application of BMD and hypersonic DNS to resolve metasurface-modified triads while maintaining a direct correlation between flow structures and spectral-space dynamics.

3.1. Fitting result of softness and validation of TDIBC

In this study, the acoustic metasurface consists of uniformly distributed spanwise rectangular cavities. The geometric parameters of these cavities were optimised (Zhao et al. Reference Zhao, Liu, Wen and Cheng2018a ) to maximise acoustic absorption performance. Regarding the parameters employed, the dimensionless half-width $b/{L_{\textit {ref}}} = 0.196$ , the unit-cell periods are $s/{L_{\textit {ref}}} = 0.52$ and depth $H/{L_{\textit {ref}}} = 1.642$ are all consistent with those in Zhao et al. (Reference Zhao, Wen, Long, Tian, Zhou and Wu2019). The geometric setting corresponds to a porosity of $\varphi = {{2b} \mathord {\left /{\vphantom {{2b} s}} \right .} s} = 0.75$ . These parameters were designed to control the second mode and have demonstrated their evident effectiveness in both LST and two-dimensional DNS (Zhao et al. Reference Zhao, Wen, Long, Tian, Zhou and Wu2019). Figure 1 compares the softness derived from a fitting procedure based on the theoretical model, which considers high-order diffraction (Zhao et al. Reference Zhao, Liu, Wen and Cheng2018a ). The fitting result was based on the method proposed by Douasbin et al. (Reference Douasbin, Scalo, Selle and Poinsot2018) using (2.7), which employs $C_0$ = −0.05 and the poles and residues shown in Appendix A. It is indicated that the fitting result agrees well with the theoretical model within a broadband frequency domain. This agreement enables the accurate simulation subject to both high- and low-frequency disturbances using the acoustic metasurface. The comparison among the cases with TDIBC, with meshed cavities and with a constant wall impedance model ( $v' = Ap'$ , and $A$ is a constant acoustic admittance), can be found in Appendix B. A good agreement is reached, validating the reliability of the TDIBC code in this study.

Figure 1. Comparison of wall softness between the reference impedance model (Zhao et al. Reference Zhao, Liu, Wen and Cheng2018a ) and the multi-pole fitting results (present).

Figure 2. ( $a$ ) Contours of optimal gain in the parameter space of the spanwise wavenumber and the angular frequency, where ${\omega _0}{L_{\textit {ref}}}/{u_\infty }$ = 0.1 and ${\beta _0}{L_{\textit {ref}}}$ = 0.8. ( $b$ ) Comparison of $N$ -factors between PSE and resolvent analysis. Here, $N$ -factor curves of PSE are shifted to be compared with resolvent analysis.

3.2. Result of stability analysis

In resolvent analysis, the interval 0–0.5 m is utilised for the optimisation problem, and the forcing is introduced at $x_0$ = 0.04 m. The $N$ -factor is defined by $N = 0.5\log ({{{E_{\textit{Chu}}}} \mathord {\left /{\vphantom {{{E_{\textit{Chu}}}} {{E_{\textit {Chu,0}}}}}} \right .} {{E_{Chu,0}}}})$ , where ${E_{\textit {Chu,0}}}$ is measured at $x$ = $x_0$ and $N$ is 0 at $x_0$ .

The gain contour for a wide range of spanwise wavenumbers and angular frequencies is shown in figure 2( $a$ ). Clearly, the low-frequency oblique wave related to the first mode ( $\omega /{\omega _0}$ , $\beta /{\beta _0}$ ) = (3, 1) and the high-frequency planar wave related to the second mode ( $\omega /{\omega _0}$ , $\beta /{\beta _0}$ ) = (10, 0) are prominent in the gain contour. In this study, $\omega /{\omega _0}$ = 10 corresponds to a dimensional frequency of 125.8 kHz, which has also been shown to be the peak frequency in the energy spectrum under the considered experiment condition (Bountin, Chimitov & Maslov Reference Bountin, Chimitov, Maslov, Novikov, Egorov, Fedorov and Utyuzhnikov2013). For brevity, the optimal response for a Fourier mode with $\omega /{\omega _0}$ = 10 and $\beta /{\beta _0}$ = 0 is referred to as mode (10, 0), and the same applies to other optimal disturbances. In the following DNS, the three-dimensional transition to turbulence will be initiated using the optimal forcings determined by resolvent analysis. As shown in figure 2( $b$ ), the optimal response captured by resolvent analysis undergoes a transient growth downstream of the forcing location, which cannot be captured in PSE when a purely modal profile is introduced at the inlet. The transient-growth region is between the forcing location and the position where the growth rate ${{\rm d}N} \mathord {\left /{\vphantom {{{\rm d}N} {{\rm d}x}}} \right .} {{\rm d}x}$ from the resolvent analysis converges with that from the PSE. For instance, the interval $x$ $ \in$ [0.04 m, 0.115 m] corresponds to the transient-growth region of Fourier mode (10, 0). The PSE results indicate that modal growth begins at approximately $x$ = 0.115 m and ends at approximately $x$ = 0.34 m for (10, 0), as shown in figure 2( $b$ ). As for the optimal oblique wave (3, 1), its $N$ -factor exceeds that of the planar wave (10, 0) downstream. This is due to a longer growth region of (3, 1), despite a lower maximum growth rate.

In the numerical simulation, both modal and non-modal growths are included. Given the broadband nature of free-stream disturbances in wind tunnels and at realistic flight conditions, both the first and second modes should be considered in the investigation of hypersonic boundary-layer instability and transition. In this study, optimal disturbances $({\omega _p}, {\beta _p})$ = (10, 0) and $({\omega _o}, {\beta _o})$ = (3, 1) were utilised to initiate the unsteady flow for cases listed in table 1. These disturbances enable complicated nonlinear mode–mode interactions and spectral broadening rapidly (Guo et al. Reference Guo, Hao and Wen2023a ). The two are also representative building blocks of the various primary instabilities with different frequencies and wavenumbers. The selection of only two of them enables the fundamental observation of the excited secondary instabilities downstream.

3.3. The TDIBC for single-frequency disturbance

In this section, a single-frequency disturbance with small amplitude is introduced in a precursor simulation to study its linear evolution individually. The verification of initialisation for the optimal forcing in the DNS is depicted in Appendix C.

Figure 3 compares the evolution of the Fourier modes (10, 0) with a solid wall and TDIBC, starting from $x$ = 0.04 or $x$ = 0.115 m. The acoustic metasurface is observed to suppress the modal growth of the second mode while enhancing the transient growth, based on the definitions of transient-growth regions in § 3.2. Note that $x$ = 0.115 m is the starting point of the modal growth of Fourier mode (10, 0), according to the PSE result in figure 2( $b$ ). When TDIBC is employed in the range of 0.04–0.2 m, the total Chu’s energy is augmented at around $x$ = 0.1m, while the wall pressure fluctuation is suppressed there in figure 3( $b$ ). This indicates that the wall pressure fluctuation may not reflect the overall impact of the metasurface. The result agrees with the numerical simulation of Wang & Zhong (Reference Wang and Zhong2012) and the wind tunnel experiment of Lukashevich et al. (Reference Lukashevich, Morozov and Shiplyuk2016), which reported that the acoustic metasurface is effective only within the unstable region of the second mode. The detailed effect of the acoustic metasurface on each primitive variable downstream and upstream of the synchronisation point of the second mode is compared in figure 4.

Figure 3. ( $a$ ) The $N$ -factor evolution and ( $b$ ) dimensionless wall pressure fluctuation for different wall boundary conditions subject to the optimal forcing of mode (10, 0).

Figure 4. Comparison of dimensionless r.m.s. magnitude between linear-stage cases using solid wall (solid line) and TDIBC (dash dot-dot line) at ( $a$ ) $x$ = 0.08 m and ( $b$ ) $x$ = 0.16 m initialised by optimal mode (10, 0). The TDIBC is applied within the range of 0.04–0.2 m.

Figure 4 compares the root-mean-square (r.m.s.) values of primitive variables between the results with a solid wall and the acoustic metasurface (realised by TDIBC). It illustrates that the r.m.s. of the wall pressure fluctuation is lower with TDIBC at $x$ = 0.08 m, which agrees with figure 3( $b$ ). However, the density and temperature are significantly enhanced near the boundary-layer edge (around $y/{L_{\textit {ref}}}$ = 2). This leads to a higher total energy of the planar wave (10, 0) in this region, namely a larger $N$ -factor shown in figure 3( $a$ ). Downstream of the neutral point of the second mode, the disturbances are completely suppressed regarding both the total energy and the wall pressure fluctuation, as shown in figure 4( $b$ ).

On the contrary, the acoustic metasurface slightly destabilises the oblique optimal disturbance (3, 1), as shown in figure 5.

Figure 5. Comparison of dimensionless pressure fluctuation at the wall for the optimal wave (3, 1) between solid-wall condition and TDIBC in the linear stage.

Based on the above finding, the starting location of TDIBC for both case 2 and case 3 was selected as $x$ = 0.115 m to avoid unnecessary amplification of unstable waves. Besides, TDIBC ends at $x$ = 0.34 m in case 3 because the second mode is stable downstream of $x$ = 0.34 m (see figure 2 $b$ ). It should be noted that no prior reports have examined the impact of the metasurface in the transient-growth stage. More research is necessary to determine how the metasurface affects the Orr/lift-up processes that underpin the transient growth. More importantly, it is not yet understood how the transition initiated by the first and second modes responds to the acoustic metasurface. In the next section, the effect of the acoustic metasurface in combination resonance (Guo et al. Reference Guo, Hao and Wen2023a ) is focused on, which requires the involvement of both first and second modes.

5.1. Evolution and contribution of the Fourier mode

The instantaneous skin friction coefficient $C_f$ is plotted in figure 10. The streamwise wavelengths of the planar wave $\lambda _{x,\,p}$ , the oblique wave $\lambda _{x,\,o}$ and large-scale $\varLambda \hbox{-}$ vortices related to the detuned mode $\lambda _{x,{\textrm { (2,2)}}}$ can be clearly identified in case 1, as shown in figure 10( $a$ ). Notably, the aligned structure of the second mode (at around $x$ = 0.2 m) loses its signature in case 2 and case 3 due to the suppression effect of the acoustic metasurface. The detuned mode (2, 2) is generated through

(5.1) \begin{align} \begin{array}{l} (3,\,1) + (3,\, - 1) \to (6,{\textrm { 0}}),\\[0.8ex] (10,{\textrm { 0}}) - (3,\, - 1) \to (7,{\textrm { 1}}),\\[0.8ex] (7,\,1) - (6,{\textrm { 0}}) \to (1,{\textrm { 1}}),\\[0.8ex] (1,\,1) + (1,\,1) \to (2,{\textrm { 2}}). \end{array} \end{align}

Figure 10. Instantaneous skin friction coefficient $C_f$ for ( $a$ ) case 1, ( $b$ ) case 2 and ( $c$ ) case 3, where $x_{\textit {onset}}$ refers to the evaluated starting location of the transition.

The generation of detuned mode (2, 2) is confirmed by wave vectors $({\alpha _r},\,\beta )$ of the considered triads (Fezer & Kloker Reference Fezer and Kloker2000), and details can be seen in Appendix D.

The detuned mode, manifested as large-scale $\varLambda \hbox{-}$ vortices (Guo et al. Reference Guo, Hao and Wen2023a ), can also be observed in case 3. This observation indicates that the flow field can be recovered as the wall boundary reverts to the solid-wall condition and that the acoustic metasurface only produces a local effect. In contrast, the detuned-mode-related $\varLambda \hbox{-}$ structures disappear in case 2, where the metasurface persists in the whole domain. In case 2, the distinct spanwise-aligned structure near $x$ = 0.45 m indicates the importance of planar wave (2, 0), as later demonstrated by figure 12( $b$ ). Furthermore, the spanwise scale approaches $0.5{\lambda _{z,\,o}}$ , which is identical to the spanwise wavelength of the detuned mode (2, 2). The enhanced pattern related to detuned modes may be responsible for the increment of skin friction in the late transitional stage with the acoustic metasurface (see case 2 in figure 7). A further decomposition of skin friction is needed to identify the dominant contributor. In the subsequent investigation, temporal and spanwise bi-Fourier analysis and energy budget analysis will be applied to highlight the key contribution to wall friction.

The Chu energy integral in the wall-normal direction includes kinetic and internal energy. Its streamwise evolution of multiple Fourier modes can comprehensively reflect the dominant modes in the transitional stage, which already include the effect of the nonlinear interaction. As depicted in figure 11, the growth of Chu’s energy of the second-mode-related mode (10, 0) is suppressed by the acoustic metasurface downstream of $x$ = 0.115 m. Meanwhile, the first mode is slightly promoted near $x$ = 0.2 m, which is consistent with previous theoretical and simulation research (Fedorov et al. Reference Fedorov2003; Wang & Zhong Reference Wang and Zhong2012). In case 1, the planar wave (10, 0) saturates at around $x$ = 0.18 m, and its amplitude begins to decrease due to energy transfer via nonlinear interaction with modes (5, 1), (10, 2) and (2, 2) or (2, 0) through the subharmonic resonance, the fundamental resonance and the combination resonance, respectively. The oblique breakdown associated with mode (3, 1), the fundamental resonance associated with mode (10, 0) and the subharmonic resonance associated with mode (10, 0) are attributable to the nonlinear interactions

(5.2) \begin{align} \begin{array}{l} (3,\,1) - (3,\, - 1) \to (0,{\textrm { 2}}),\\[0.8ex] (10,{\textrm { 0}}) + (0,{\textrm { 2}}) \to (10,{\textrm { 2}}),\\[0.8ex] (10,{\textrm { 0}}) - (5,{\textrm { 1}}) \to (5,{ - 1}), \end{array} \end{align}

respectively. The aforementioned modes other than modes $(3, \pm 1)$ and (10, 0), known as secondary instabilities, undergo a rapid growth process and make the flow develop into the late stage of transition and breakdown. The specific nonlinear mode–mode interaction mechanism was systematically studied and reported by Guo et al. (Reference Guo, Hao and Wen2023a ). In this study, the effect of the acoustic metasurface is focused on. In figure 11( $a$ ), one can see that the second-mode-related mode (10, 0) in case 2 (dashed line) grows gradually to nearly the same amplitude as that in case 1 (solid line), but saturates further downstream. Meanwhile, the secondary growth of the primary oblique wave (3, 1) shows a lower growth rate (i.e. slope) in case 2 at around $x$ = 0.28 m than in case 1. This result is somewhat consistent with the study of Chen et al. (Reference Chen, Zhu and Lee2017), where the amplitude of the second mode saturates and then the ‘phase lock’ sets in. This phase lock supports a secondary growth stage of the primary low-frequency mode. Meanwhile, the second mode decays due to energy transfer to the rapidly growth secondary modes. Regarding the metasurface effect in this paper, it suppresses the first-mode-related (3, 1) when the nonlinearity is considered. This finding differs from the linear case in figure 5, which is probably due to less energy being transferred from the attenuated mode (10, 0) with the metasurface. Another interesting finding is that, the dominant detuned modes (2, 0) are apparently more energetic in case 2 and case 3, compared with the solid-wall case. In addition, the low-frequency modes generated in the later transitional stage, such as (1, 0), are also more evident than their counterparts in the solid-wall case.

Figure 11. Comparison of streamwise development of Chu’s energy: ( $a$ ) and ( $b$ ) for case 1 and case 2, and ( $c$ ) and ( $d$ ) for case 1 and case 3. The modes in case 1 are represented by solid lines with filled symbols, and modes in case 2 and case 3 are represented by dotted lines with open symbols.

Figure 12. Comparison of maximum absolute modal contribution to the instantaneous skin friction coefficient $C_f$ : ( $a$ ) and ( $b$ ) for case 1 and case 2, and ( $c$ ) and ( $d$ ) for case 1 and case 3. The modes in case 1 are represented by solid lines with filled symbols, and modes in case 2 and case 3 are represented by dotted lines with open symbols.

We intend to identify the dominant Fourier mode contributing to the increment of higher skin friction in case 2. To this end, the maximum absolute contribution from mode ( $m$ , $n$ ) to the instantaneous skin friction is defined by

(5.3) \begin{align} \Delta {C_{f,(m,n)}}(x) = \mathop {\max }\limits _{(z,t)} \left | {{C_{f,(m,n),\textit {unsteady}}} - {C_{f,\textit {laminar}}}} \right | \!, \end{align}

where $C_{f,(m,n),\textit {unsteady}}$ and $C_{f,\textit {laminar}}$ are the instantaneous skin friction induced by the laminar flow plus mode ( $m$ , $ \pm n$ ) and by the laminar flow alone, respectively. This indicator is defined to highlight the Fourier modes which are the source contributors to the skin friction in figure 10. Figure 12 compares the maximum absolute modal contribution to the instantaneous skin friction coefficient $C_f$ between cases with (case 2 and case 3) and without (case 1) the acoustic metasurface. Apparently, the second mode (10, 0) is suppressed in both case 2 and case 3 with the acoustic metasurface, which agrees with the disappearence of the aligned structure in the vortex visualisation (shown in figure 6) and the instantaneous wall friction (shown in figure 10). As displayed in figure 12( $b$ ), the detuned modes (2, 2) and (2, 0) demonstrate a considerably higher contribution, indicating a significantly strengthened combination resonance. In addition, the subharmonic-related mode (5, 1) and oblique mode (3, 1) in case 2 also show a more pronounced contribution than those in case 1. This result indicates the promotion of these two nonlinear interaction mechanisms (combination and subharmonic resonance) using the acoustic metasurface. Note that there are some differences between the wall friction contributions and the integrated Chu energy, which resembles the linear case in § 3.3. For example, the Chu energies of modes (2, 2) and (2, 0) are decreased and increased by the metasurface in figure 11, respectively. However, the contributions to $C_f$ of both modes are notably augmented by the metasurface. This implies that the wall quantity and the integrated energy across the boundary layer may have different performances. To quantitatively evaluate the effect of the metasurface, the difference of $\Delta C_f$ is calculated between case 1 and case 2. The result serves to identify the main contributor to the strengthened skin friction in the late transitional region in case 2. The overall integral effect is calculated by

(5.4) \begin{align} {\delta _{(m,\:n)}} = \int _{{x_1}}^{{x_2}} {\left ( {\Delta {C_{f,\:(m,n),\:\textit {Case 2}}} - \Delta {C_{f,\:(m,n),\:\textit {Case 1}}}} \right )\,{\rm d}x}, \end{align}

where ${x_1}$ = 0.34 m and ${x_2}$ = 0.6 m. The comparison of $\delta _{(m,\:n)}$ , normalised by $\delta _{(0,\:0)}$ , among different Fourier modes is presented in figure 13. The results indicate that the dominant contributor is the detuned planar mode (2, 0), which aligns with the predominant planar structure observed at approximately $x$ = 0.45 m in figure 10( $b$ ). The above analysis answers the question raised by the second research objective in the Introduction.

Figure 13. Comparison of $\delta _{(m,\; n)}$ among different Fourier modes. The value $\delta _{(0, \:0)}$ is utilised for normalisation.

Regarding the question of why transition onset is delayed, the increasing contribution of the MFD (0, 0) should be discussed. In the range $x \in [0.25, 0.32]$ m where the transition begins in figure 7, the presence of the metasurface leads to the postponed growing contribution of (0, 0) in figure 12( $a$ ). In this range, only the streak mode (0, 2) is notable in terms of the magnitude. As a result, the delayed transition onset can be attributed to the postponed growth of mode (0, 2). This finding is interesting because the manifested control effect by the metasurface does not lie in the second mode directly. Given the nonlinearity, the Fourier mode (0, 2) is the direct contributor to the delayed transition onset instead. Upstream of $x$ = 0.3 m, the contributions of other secondary modes are not evident in figure 12( $a$ ). As a result, the delayed (0, 2) is mainly caused by the suppressed oblique-wave resonance (see the first interaction equation in (5.2)) via the relatively large-amplitude primary instability wave (3, 1). As interpreted in the discussion of figure 11, mode (3, 1) is stabilised by the metasurface in the nonlinear stage, which is probably achieved by the less energy received from the suppressed Mack second mode. Therefore, the second mode plays an indirect role in the transition delay mechanism of the acoustic metasurface.

5.2. Mode–mode interaction

Through complementary spectral analyses, the bi-Fourier decomposition quantifies disturbance amplitude development, while the bispectral analysis elucidates the underlying energy-transfer pathways governing these evolutionary patterns. To clarify the energy-transfer mechanism, BMD was performed on the mid- $xy$ plane of the computational domain, where the dominant triadic interactions are clearly discernible through the bispectrum value distribution. As illustrated in figure 14, the bispectrum maps at various streamwise sections reveal distinct nonlinear interaction patterns, with the magnitude of bispectrum values $\phi$ directly correlating with the strength of triadic coupling. The upper half- and lower half-planes of each map correspond to sum and difference triadic interactions, respectively. The coordinates ( $k$ , $l$ ) and ( $m$ , - $l$ ) represent the sum interaction ( ${\omega _k} + {\omega _l} \to {\omega _m}$ ) and the difference interaction ( ${\omega _m} - {\omega _l} \to {\omega _k}$ ). Here, the subscripts $k$ , $l$ , $m$ denote non-dimensional angular frequencies (e.g. $\omega _{10}$ represents Fourier modes of $\omega$ = 10, equivalent to the physical frequency $f$ = 125.8 kHz) and satisfy the triad resonant condition $k + l = m$ . Diagonal dashed lines demarcate triadic interactions generating constant frequencies, as annotated by their respective angular frequencies $\omega$ . These dashed lines in the figures cover all the triad interactions that generate the mode with the corresponding angular frequency $\omega$ . The initial stage exhibits prominent interactions between the high-frequency component ( $\omega _{10}$ ) and low-frequency mode ( $\omega _{3}$ ), generating $\omega _{7}$ and $\omega _{13}$ modes through difference ( ${\omega _{10}} - {\omega _3} \to {\omega _7}$ ) and sum ( ${\omega _{10}} + {\omega _3} \to {\omega _{13}}$ ) triadic interactions, respectively, as displayed in figure 14( $a{-}b$ ). Concurrently, sum-self-interaction of $\omega _{3}$ produces its harmonic $\omega _{6}$ , while subsequent difference interaction between $\omega _{7}$ and $\omega _{3}$ yields $\omega _{4}$ . Downstream evolution demonstrates progressive attenuation of high-frequency components accompanied by amplification of low-frequency modes, particularly the detuned $\omega _{2}$ and $\omega _{1}$ modes generated through ${\omega _{3}} - {\omega _1} \to {\omega _2}$ and ${\omega _{3}} - {\omega _2} \to {\omega _1}$ interactions (see figure 14 $f$ ). The late transitional region ( $x \in [0.55, 0.575]$ ) exhibits emergent near-zero-frequency components with the absence of dominant spectral peaks, indicative of fully developed turbulent flows. This energy-transfer pathway, where triadic interactions mediate progressive energy transfer from high- to low-frequency modes culminating in near-zero-frequency components, contrasts with previous hypersonic boundary-layer studies employing bicoherence analysis. While Craig et al. (Reference Craig, Humble, Hofferth and Saric2019) and Zhang & Shi (Reference Zhang and Shi2022) reported dominant self-sum interactions of the second mode ( $\omega _{10}$ ) and its harmonic in quiet wind tunnel conditions, the present results highlight stronger detuned mechanisms involving first-mode interactions. This discrepancy likely stems from the relative suppression of second-mode self-sum interactions under the current configuration and flow condition, combined with enhanced first-mode activity that is typically attenuated in the quiet wind tunnel experimental environments.

Figure 14. Magnitude of mode bispectrum of case 1 for various streamwise regions ( $a$ ) 0.2–0.225, ( $b$ ) 0.25–0.275, ( $c$ ) 0.3–0.325, ( $d$ ) 0.35–0.375, ( $e$ ) 0.4–0.425, ( $f$ ) 0.45–0.475, ( $g$ ) 0.5–0.525 and ( $h$ ) 0.55–0.575 m at $z$ = $z_0$ /2, where $z_0$ is the spanwise width of the computational domain.

To confirm the spanwise wavenumber for the dominant triadic interactions, the spatial structure of bispectral modes was revealed using BMD. Figure 15 presents the $xz$ -plane bispectral modes characterised by density at three characteristic wall-normal positions: the wall surface ( $y$ = 0 mm), an intermediate layer ( $y$ = 1.6 mm) and the boundary-layer edge vicinity ( $y$ = 3.2 mm). The results reveal that the triadic interaction associated with detuned modes (corresponding to $\omega _{1}$ and $\omega _{2}$ ) exhibits high energy density near the wall, whereas the first mode ( $\omega _{3}$ ) and its harmonic ( $\omega _{6}$ ) are more energetic in the off-wall region (based on the bispectrum map at three characteristic wall-normal positions, not shown here). This spatial distribution is consistent with both the vortical structure characteristics of the first instability mode near the boundary-layer edge and the dominant influence of detuned modes on skin friction generation during late transitional stages, as documented in Guo et al. (Reference Guo, Hao and Wen2023a ). Notably, the spanwise bispectral mode distribution serves as an effective indicator of the dominant spanwise wavenumber for specific frequencies. For instance, modes (1, 1), (3, 1), (6, 0) and (7, 1) are identified as primary contributors to $\omega _{1}$ , $\omega _{3}$ , $\omega _{6}$ and $\omega _{7}$ , respectively. Modes (2, 2) and (2, 0) jointly govern $\omega _{2}$ formation, with (2, 0) exhibiting greater dominance due to the weak spanwise periodicity (see figure 15 $j$ ). These detuned modes are further responsible for the large-scale Lambda structures observed near $x$ = 0.4 m in figure 15. In summary, the BMD successfully resolves the dominant spanwise features of the transitional flow dynamics, involving both primary (first and second) instability modes as well as secondary modes. The triadic interaction mechanisms are quantitatively elucidated during the transition process.

Figure 15. Bispectral mode (real part of density component) marked with the associated sum and difference triadic interactions at $y$ = 3.2 ( $a{-}c$ ), $y$ = 1.6 ( $d{-}h$ ) and $y$ = 0 mm ( $i{-}k$ ).

Figure 16 presents the bispectrum distribution at different streamwise locations for case 2. To highlight the influence of the acoustic metasurface on nonlinear mode–mode interactions, the colour bar levels are maintained consistent with the solid-wall case in figure 14 for a direct comparison. The introduction of the acoustic metasurface (case 2) leads to a pronounced suppression of second-mode-dominated triadic interactions in the region $x \in [0.25\,\mathrm{m}, 0.275\,\mathrm{m}]$ (see figure 16 $b$ , compared with figure 14 $b$ ). This suppression subsequently weakens the triadic interaction between $\omega _{6}$ and $\omega _{3}$ in $x \in [0.3\,\mathrm{m}, 0.325\,\mathrm{m}]$ (see figure 16 $c$ , compared with figure 14 $c$ ), owing to reduced energy transfer from $\omega _{10}$ . This observation also echoes with the deduction in § 5.1. As the generation triad for $\omega _3$ is suppressed, the resulting stationary streak mode is less energetic and the MFD growth and transition onset in figure 12 are thus delayed by the metasurface. Further downstream, the amplified presence of $\omega _{1}$ and $\omega _{2}$ (see figure 16 $e{-}g$ , compared with figure 14 $e{-}g$ ) indicates stronger nonlinear interactions involving detuned modes with the presence of the metasurface. In other words, the metasurface effect on the low-frequency dynamics is evidently different, where the low-frequency mode interactions are enhanced in the controlled case. This finding is consistent with the elevated skin friction observed for case 2 beyond $x = 0.34\,\mathrm{m}$ (figure 7). The enhanced triadic interactions generating the detuned modes result in an increased contribution to the skin friction coefficient, as shown in figure 12. This effect may originate from the shear-induced dissipation observed in figure 9 from the perspective of mean-flow statistics. Notably, the underlying triadic interaction mechanism remains unaltered, as low-frequency detuned modes ( $\omega _{2}$ and $\omega _{1}$ ) continue to emerge in later transitional stages. Interestingly, residual light points persist in the triadic maps for case 2 (figure 16 $h$ ), suggesting that the flow in the controlled case may not yet be fully developed by the end of the computational domain. Further analysis of the flow state will be provided in subsequent sections.

Figure 16. Magnitude) of mode bispectrum of case 2 for various streamwise regions ( $a$ ) 0.2–0.225, ( $b$ ) 0.25–0.275, ( $c$ ) 0.3–0.325, ( $d$ ) 0.35–0.375, ( $e$ ) 0.4–0.425, ( $f$ ) 0.45–0.475, ( $g$ ) 0.5–0.525 and ( $h$ ) 0.55–0.575 m at $z$ = $z_0$ /2, where $z_0$ is the spanwise width of the computational domain.

5.3. Paths to turbulent flow

In this section, the spanwise- and time-averaged flow field is focused on and analysed. To evaluate the flow state and its progression, the van-Driest-transformed velocity (Van Reference Van and Edward1951) $U_{{\textit{VD}}}^ +$ profiles for three cases at various streamwise stations are plotted in figure 17. The influence of mean density variations due to the effect of compressibility can be removed by van Driest transformation for the mean velocity

(5.5) \begin{align} U_{{\textit{VD}}}^ + = \int _0^{{{\bar u}^ + }} {\sqrt {\frac {{\bar \rho }}{{{{\bar \rho }_w}}}} } {\rm d}{\bar u^ + }(y), \end{align}

where the friction velocity $\bar u_\tau$ is utilised for normalisation

(5.6) \begin{align} {\bar u^ + } = \frac {{\bar u}}{{{{\bar u}_\tau }}}, \: \,{\bar u_\tau } = \sqrt {\frac {{{{\bar \tau }_w}}}{{{{\bar \rho }_w}}}}. \end{align}

Figure 17. Van Driest velocity profiles for case 1, case 2 and case 3 at ( $a$ ) $x$ = 0.2, ( $b$ ) $x$ = 0.3, ( $c$ ) $x$ = 0.4, ( $d$ ) $x$ = 0.5, ( $e$ ) $x$ = 0.55, ( $f$ ) $x$ = 0.6 m.

The viscous sublayer law ( ${u^ + } = {y^ + }$ ) and the log-layer law ( ${u^ + } = 5.8 + {1 \mathord {\left /{\vphantom {1 \kappa }} \right .} \kappa }\ln {y^ + }$ , $\kappa$ = 0.4) are given as references, as black dashed lines. The intercept 5.8 is larger than the incompressible counterpart 5.0, which has been reported by Guo et al. (Reference Guo, Shi, Gao, Jiang, Lee and Wen2022a ). At $x$ = 0.2 m, the result of the three cases agrees well with the reference, showing a consistent and well-resolved viscous sublayer. Downstream $x$ = 0.3 m, case 1 shows its initial establishment of the log layer, while in case 2 and case 3, velocity profiles just begin to deviate from sublayer reference at around $y^ +$ = 10. Starting from $x$ = 0.55 m, the velocity profile gradually approaches the reference log-layer law, and the transformed velocity profile is nearly the same for case 1 and case 3 at $x$ = 0.6 m. The observations indicate that fully developed turbulence is established, even if the slope shows a little difference due to the high Mach number. A similar difference was also reported in supersonic/hypersonic simulation by Zhou et al. (Reference Zhou, Liu, Lu and Yan2023) and Guo et al. (Reference Guo, Shi, Gao, Jiang, Lee and Wen2022a ). As for case 2, the near-wall slope ( ${{\rm d}U_{{\textit{VD}}}^ + } \mathord {\left /{\vphantom {{dU_{{\textit{VD}}}^ + } {d{y^ + }}}} \right .} {{\rm d}{y^ + }}$ ) remains lower due to the acoustic metasurface utilised, compared with the solid-wall case downstream $x$ = 0.3 m (see figure 17 $c{-}f$ ). This is qualitatively consistent with the permeable-wall effect on the supersonic boundary layer studied by Chen & Scalo (Reference Chen and Scalo2021a , Reference Chen and Scalo2021b ).

Owing to the augmented near-wall shear-induced dissipation in the late transitional region for case 2, the van-Driest-transformed velocity profiles deviate from the viscous law near the wall. To evaluate the flow state, the PSD of wall pressure fluctuations at the downstream location is analysed. As illustrated in figure 18, a good agreement is reached with the reference-scale law for case 1 at $x$ = 0.6 m, indicating the establishment of a fully developed turbulent state. In contrast, the comparison with the $\omega ^{ - 5/3}$ law in the inertial subrange for case 2 shows less consistency compared with case 1. Moreover, the peak frequency ( $f$ = 25.2 kHz) corresponding to detuned modes (2, 0) and (2, 2) is observed at $x$ = 0.6 m for case 2. This indicates that spectral broadening is incomplete at this location, and fully developed turbulence has not yet been established.

Figure 18. Power spectral density (PSD) of wall pressure fluctuation at $x$ = 0.45 m for ( $a$ ) case 1 and ( $b$ ) case 2, and at $x$ = 0.6 m for ( $c$ ) case 1 and ( $d$ ) case 2, in comparison with the law of $\omega ^{ - 5/3}$ in the inertial subrange and the law of $\omega ^{ - 7}$ in the dissipation scale. The green dash-dot-dot line marks the frequency corresponding to detuned modes (2, 0) and (2, 2).

In summary, the transition onset is delayed by the acoustic metasurface, and fully developed turbulence, as observed in the solid-wall case, does not appear to have been established within the computational domain using the acoustic metasurface in case 2.