1. Introduction
Shock-wave–boundary-layer interactions (SWBLIs) are commonly encountered in practical applications such as transonic wings, inlets of supersonic engines and compressor blades in turbomachinery (Sabnis & Babinsky Reference Sabnis and Babinsky2023). These interactions pose severe challenges for high-speed vehicles, leading to issues such as aeroheating overload, structural fatigue and engine inlet unstart. Moreover, SWBLIs significantly impact flow instability, flow dynamics and transition processes. Consequently, SWBLIs are of considerable value in both applications of engineering and fundamentals of fluid mechanics (Dolling Reference Dolling2001; Datta Reference Datta2015).
Since first documented in experimental works in the 1940s (Liepmann Reference Liepmann1946; Ackeret, Feldmann & Rott Reference Ackeret, Feldmann and Rott1947), SWBLIs have received sustained attention and research (Gadd, Holder & Regan Reference Gadd, Holder and Regan1954; Reyhner & Flügge-Lotz Reference Reyhner and Flügge-Lotz1968; Adamson Jr & Messiter Reference Adamson and Messiter1980; Degrez, Boccadoro & Wendt Reference Degrez, Boccadoro and Wendt1987; Stanewsky Reference Stanewsky1988; Katzer Reference Katzer1989; Settles & Dodson Reference Settles and Dodson1991). Most of these studies have been limited to examining the effects of Mach number, Reynolds number and shock intensity on steady aerodynamic quantities such as surface pressure, skin friction and aerodynamic heating (Dolling Reference Dolling2001). Since the 1990s, the importance of investigating the instabilities in laminar high-speed shear layers and boundary layers was recognized (Fedorov & Khokhlov Reference Fedorov and Khokhlov1991; Sandham & Reynolds Reference Sandham and Reynolds1991; Balakumar & Malik Reference Balakumar and Malik1992; Kachanov Reference Kachanov1994; Hanifi, Schmid & Henningson Reference Hanifi, Schmid and Henningson1996), leading to a major focus on SWBLI-induced instabilities and transition.
Classification of the instabilities in SWBLIs.
Flow instabilities influenced by SWBLIs can be categorized into modal and non-modal types, as illustrated in figure 1. Modal instabilities correspond to discrete eigenmodes of the linearized Navier–Stokes (LNS) operator, typically characterized by exponential growth or decay in time or space. In contrast, non-modal instability arises from the non-normality of the linearized operator, which enables transient amplification of disturbances without unstable eigenvalues. Rather than growing exponentially, non-modal mechanisms amplify disturbances algebraically, through transfer of energy from the two cross-streams to the streamwise velocity component, a process known as the lift-up effect (Ellingsen & Palm Reference Ellingsen and Palm1975; Landahl Reference Landahl1980). Both the modal and non-modal instabilities can trigger the transition, typically through the breakdown of streaks, which is observed experimentally (Sandham et al. Reference Sandham, Schülein, Wagner, Willems and Steelant2014; Giepman, Schrijer & Van Oudheusden Reference Giepman, Schrijer and Van Oudheusden2015; Willems, Gülhan & Steelant Reference Willems, Gülhan and Steelant2015; Currao et al. Reference Currao, Choudhury, Gai, Neely and Buttsworth2020; Butler & Laurence Reference Butler and Laurence2022; Benitez et al. Reference Benitez, Borg, Scholten, Paredes, McDaniel and Jewell2023, Reference Benitez2025; Mahalingesh, Piponniau & Dupont Reference Mahalingesh, Piponniau and Dupont2023) and numerically (Novikov Reference Novikov2017; Currao et al. Reference Currao, Choudhury, Gai, Neely and Buttsworth2020; Fu et al. Reference Fu, Karp, Bose, Moin and Urzay2021; Lugrin et al. Reference Lugrin, Beneddine, Leclercq, Garnier and Bur2021; Dwivedi, Sidharth & Jovanović Reference Dwivedi, Sidharth and Jovanović2022; Cao et al. Reference Cao, Hao, Klioutchnikov, Wen, Olivier and Heufer2022).
Modal instabilities can be further divided into local (convective) instability and global instability. Local instability, often referred to as convective instability, is characterized by spatial amplification of perturbations that are convected downstream. Mack’s modes (Mack Reference Mack1984), such as the first and second modes, are representative examples of convective instability. Linear stability theory (LST) (Malik Reference Malik1989, Reference Malik1997; Reed, Saric & Arnal Reference Reed, Saric and Arnal1996) is an effective method for capturing Mack’s modes, especially when the flow is weakly non-parallel. In the eigenvalue spectrum obtained from LST, three primary branches are typically identified: the fast acoustic branch; the slow acoustic branch; the entropy–vorticity branch (Fedorov Reference Fedorov2011). The fast acoustic branch corresponds to acoustic disturbances with phase speeds higher than the free stream velocity and is mainly associated with pressure fluctuations propagating outside the boundary layer. The slow acoustic branch represents acoustic waves with lower phase speeds, often localized near the boundary layer. In contrast, the entropy–vorticity branch is characterized by disturbances dominated by entropy and vorticity fluctuations, convecting with the mean flow inside the boundary layer. Combined with the development of direct numerical simulation (DNS) and large-eddy simulation (LES), the destabilization of Mack’s modes and their effects on the transition process have been explored in detail. For example, Fasel et al. (Reference Fasel, Thumm, Bestek, Kral and Zang1993) introduced a pair of oblique first modes in a Mach 3 flat plate boundary layer, demonstrating that oblique breakdown could lead to a fully turbulent boundary layer. The entire process, from linear, and nonlinear growth, to eventual breakdown, was investigated numerically using LST, parabolized stability equations (PSE) and DNS. Chang & Malik (Reference Chang and Malik1994) also simulated the oblique breakdown and further discussed the secondary instability. In the context of SWBLIs, Adams (Reference Adams2001) showed that the second mode is stable in the separated boundary layer along a compression corner at Mach 5. Additionally, Balakumar, Zhao & Atkins (Reference Balakumar, Zhao and Atkins2005) performed LST and DNS on a
$5.5^{\circ }$
compression corner under Hyper-X wind-tunnel condition, revealing that higher Mack’s modes appear in the separation region but cannot sustain growth across the entire separated region. Yao et al. (Reference Yao, Krishnan, Sandham and Roberts2007) investigated the impact of Mach number on the most unstable mode in SWBLIs using LST and PSE.
Görtler instability is another type of convective-type instability. Its growth is caused by the concave surface in boundary layer flow, manifesting as streamwise-oriented, counter-rotating vortices (Ginoux Reference Ginoux1971; Floryan Reference Floryan1991; Saric Reference Saric1994). These vortices can lead to the formation of low- and high-speed streaks. In the flat plate boundary layer without curvature, streaks are typically associated with component-type instability (Hanifi et al. Reference Hanifi, Schmid and Henningson1996; Bugeat et al. Reference Bugeat, Chassaing, Robinet and Sagaut2019), originating from transient growth due to the lift-up mechanism. In SWBLIs, streaks have been widely observed downstream of reattachment experimentally (Heffner, Chpoun & Lengrand Reference Heffner, Chpoun and Lengrand1993; de la Chevalerie et al. Reference de la Chevalerie, Fonteneau, De Luca and Cardone1997; Bleilebens & Olivier Reference Bleilebens and Olivier2006; Roghelia et al. Reference Roghelia, Olivier, Egorov and Chuvakhov2017; Lugrin et al. Reference Lugrin, Nicolas, Severac, Tobeli, Beneddine, Garnier, Esquieu and Bur2022), accompanied with significantly increased peak heating and spanwise variations. The amplification of streaks in SWBLIs can be attributed to various mechanisms, including baroclinic effects (Dwivedi et al. Reference Dwivedi, Sidharth, Nichols, Candler and Jovanović2019), nonlinear interactions (Lugrin et al. Reference Lugrin, Beneddine, Leclercq, Garnier and Bur2021) and the saturation of global instabilities (Cao et al. Reference Cao, Hao, Klioutchnikov, Olivier and Wen2021). Additionally, Görtler instability induced by the concave curvature formed at the separation and reattachment points is another significant contributor (Navarro-Martinez & Tutty Reference Navarro-Martinez and Tutty2005; Chuvakhov et al. Reference Chuvakhov, Borovoy, Egorov, Radchenko, Olivier and Roghelia2017; Roghelia et al. Reference Roghelia, Olivier, Egorov and Chuvakhov2017).
Resolvent analysis (Ehrenstein & Gallaire Reference Ehrenstein and Gallaire2005; Alizard & Robinet Reference Alizard and Robinet2007; Monokrousos et al. Reference Monokrousos, Åkervik, Brandt and Henningson2010), more plainly called input–output analysis (Dwivedi et al. Reference Dwivedi, Sidharth, Nichols, Candler and Jovanović2019), takes the component-type non-normality into account, making it capable of investigating both convective-type and component-type instabilities. Moreover, it is a fully non-parallel method that does not require any assumption of the base flow. Resolvent analysis has been applied in various studies to investigate flow instabilities and related mechanisms. Specifically, Bugeat et al. (Reference Bugeat, Chassaing, Robinet and Sagaut2019) utilized it to capture streaks, the first mode and the second mode in a boundary layer at Mach 4.5. Dwivedi et al. (Reference Dwivedi, Sidharth, Nichols, Candler and Jovanović2019) focused on streaks exclusively and suggested that the streamwise deceleration and baroclinic effect of the separation bubble are responsible for the amplification of streaks. Recently, Hao et al. (Reference Hao, Cao, Guo and Wen2023) conducted resolvent analysis on a Mach 7.7 compression corner and concluded that the growth of streaks should be attributed to the Görtler instability when the flow is globally stable.
Despite these studies, the amplification mechanism of streaks by SWBLIs remains inconclusive, and limited research has been carried out on the amplification mechanism of Mack’s modes, particularly the first mode. Dwivedi et al. (Reference Dwivedi, Sidharth and Jovanović2022) conducted resolvent analysis on a double-wedge configuration and demonstrated that the amplification of oblique waves arises from the shear layer due to curvature. Benitez et al. (Reference Benitez, Borg, Scholten, Paredes, McDaniel and Jewell2023) conducted experiments on a Mach 6 cone–cylinder–flare and identified the second mode and a ‘shear-layer’ instability in the same frequency band as the first mode. For the same configuration and Mach number, Caillaud et al. (Reference Caillaud, Lugrin, Esquieu and Content2023) presented resolvent analysis results of streaks, the first mode and the second mode. A further investigation by Caillaud et al. (Reference Caillaud, Sandham, Touber and Sandberg2025) used LST and inviscid energy budgets to highlight the characteristics in the low- and high-frequency responses, corresponding to the first and second modes in nature, respectively. However, the mechanisms underlying how SWBLIs amplify flow instabilities remain poorly understood. Although Hao et al. (Reference Hao, Cao, Guo and Wen2023) provided insights into the amplification streaks and the second mode at Mach 7.7, the first mode was absent due to the relatively cold wall of
$T_w/T_0=0.18$
. Consequently, the distinction between the first mode and the so-called ‘shear-layer’ instability is blurred, leaving its amplification mechanisms undetermined. The above review highlights one of the main focuses of this paper: revealing the amplification mechanisms of these instabilities in SWBLIs, including streaks, the first mode and the second mode. A key aspect is the first mode, where its amplification mechanism and underlying stability explanation will be thoroughly examined using a combination of resolvent analysis and LST.
Referring back to the categories in figure 1, as the strength of SWBLI is sufficiently strong, a larger separation bubble can support global instability. This is an intrinsic instability of the separation bubble that does not require any external perturbation. In the incipient stage of global instability, only a stationary unstable mode is present. With a further increase in SWBLI strength, additional oscillatory modes appear. Robinet (Reference Robinet2007) identified the existence of a stationary three-dimensional (3-D) globally unstable mode in an incident shock flow at Mach 2.15 through global stability analysis (GSA). Similarly, Hildebrand et al. (Reference Hildebrand, Dwivedi, Nichols, Jovanović and Candler2018) conducted GSA on an oblique SWBLI at Mach 5.92, finding this stationary global mode and attributing the downstream streaks to the global instability. Hao et al. (Reference Hao, Cao, Wen and Olivier2021) investigated the characteristics of global instability at various angles and temperatures and established a critical criterion for compression ramp. Song & Hao (Reference Song and Hao2023) performed GSA and DNS on an oblique SWBLI, establishing a similar critical criterion for shock incidence and observing unsteadiness originating from a stationary global mode. Li & Hao (Reference Li and Hao2023) further extended the critical criterion to a hollow cylinder–flare. Most of the above studies on global instability focus on a single interaction type, despite the similarities between compression corner and incident shock flows. This paper aims to compare the global characteristics of compression corner and incident shock flows, highlighting their commonalities and differences.
Based on the above review, this paper seeks to depict a comprehensive picture of instabilities in laminar SWBLIs, by means of GSA, resolvent analysis and LST. Given that the compression corner and incident shock flows are often considered equivalent, the analysis is conducted for both configurations to identify their commonalities and differences. The flow conditions are based on the experiment conducted by Willems et al. (Reference Willems, Gülhan and Steelant2015), with a Mach number of 6 to cover the spectrum of instabilities, including streaks, the first mode and the second mode. The remaining content of this paper is organized as follows: § 2 introduces the numerical methods, including the base flow solver and stability analysis tools; § 3 provides the computational details; § 4 presents results in detail, corresponding to GSA and resolvent analysis; the main conclusions are summarized in § 5.
2. Numerical methods
2.1. Base flow solver
The governing equations for base flow simulations are the 3-D compressible Navier–Stokes equations, expressed in tensor form as
\begin{equation} \frac {\partial {\boldsymbol{U}}}{\partial t} + \frac {\partial {\boldsymbol{F}}_{\kern-2pt j}}{\partial x_j} = \frac {\partial {\boldsymbol{F}}_{v,j}}{\partial x_j}, \end{equation}
where
\begin{equation} \boldsymbol{U}=\left (\!\!\begin{array}{c} \rho \\ \rho u_{i}\\ \rho e \end{array}\!\!\right)\!, \quad \boldsymbol{F}_{\kern-2pt j}=\left (\!\!\begin{array}{c} \rho u_{j} \\[3pt] \rho u_{i} u_{j}+p \delta _{\textit{ij}} \\[3pt] (\rho e+p) u_{j} \end{array}\!\!\right )\!, \quad \boldsymbol{F}_{v, j}=\left (\!\!\begin{array}{c} 0 \\ \tau _{\textit{ij}} \\ \tau _{\textit{ij}} u_{i}+q_{j} \end{array}\!\!\right ) \end{equation}
are the vectors of conservative variables, inviscid fluxes and viscous fluxes, respectively. More details about the governing equations are described in Song & Hao (Reference Song and Hao2023). The base flow simulations are performed using an in-house, multiblock parallel finite-volume solver called PHAROS (Hao, Wang & Lee Reference Hao, Wang and Lee2016; Hao & Wen Reference Hao and Wen2020). The inviscid fluxes are calculated using the Harten-Lax-van Leer-Contact scheme (Harten, Lax & Leer Reference Harten, Lax and Leer1983; Toro, Spruce & Speares Reference Toro, Spruce and Speares1994), and the viscous fluxes are discretized by the second-order central scheme. For time iteration, an implicit line relaxation method (Wright, Candler & Bose Reference Wright, Candler and Bose1998) is employed to accelerate calculations.
2.2. Global stability analysis
In stability analysis, the vector of conservative variables
$\boldsymbol{U}$
is decomposed into a two-dimensional (2-D) steady solution
$\boldsymbol{U}_{2-D}$
and a 3-D small amplitude perturbation
$\boldsymbol{U}^{\boldsymbol{\prime}}$
as
\begin{equation} {\boldsymbol{U}}(x,y,z,t)={\boldsymbol{U}}_{2-D}(x,y,z,t)+{\boldsymbol{U}^{\boldsymbol{\prime}}}(x,y,z,t). \end{equation}
Substituting (2.3) into (2.1) and neglecting the high-order terms lead to the LNS equations,
\begin{equation} \frac {\partial {\boldsymbol{U}^{\boldsymbol{\prime}}}}{\partial t} + \frac {\partial {\boldsymbol{F'}}_{\kern-2pt j}}{\partial x_j} = \frac {\partial {\boldsymbol{F'}}_{\kern-1pt v,j}}{\partial x_j}, \end{equation}
in the operator form as
\begin{equation} \frac {\partial {\boldsymbol{U}^{\boldsymbol{\prime}}}}{\partial t}=\boldsymbol{A}{\boldsymbol{U}^{\boldsymbol{\prime}}}. \end{equation}
Here
$\boldsymbol{A}$
is the Jacobian matrix determined by the base flow. When the perturbation is assumed to be periodic in the spanwise direction,
$\boldsymbol{U'}$
can be expressed in the modal form as
\begin{equation} {\boldsymbol{U}^{\boldsymbol{\prime}}}(x,y,z,t)={\boldsymbol{\hat {U}}}(x,y)\exp \left [\mathrm{i}{\beta }z-\mathrm{i}(\omega _r+\mathrm{i}\omega _i)t\right ]. \end{equation}
Here
$\boldsymbol{\hat {U}}(x,y)$
is the 2-D eigenfunction,
$\beta$
is the spanwise wavenumber,
$\omega _{r}$
and
$\omega _i$
are the angular frequency and the growth rate, respectively. Substituting (2.6) into (2.5) turns the global stability problem into an eigenvalue problem as
\begin{equation} \boldsymbol{A}(\beta )\boldsymbol{\hat {U}}=-\mathrm{i}(\omega _r+\mathrm{i}\omega _i)\boldsymbol{\hat {U}}. \end{equation}
The Jacobian matrix
$\boldsymbol{A}$
is constructed from LNS using a second-order scheme. Specifically, the inviscid Jacobians are computed by the modified Steger–Warming scheme (MacCormack Reference MacCormack2014) near discontinuity and a central scheme in smooth regions. The viscous Jacobians are evaluated by the second-order central difference scheme. The detailed discretization schemes are described in Hao et al. (Reference Hao, Cao, Wen and Olivier2021). The eigenvalues of
$\boldsymbol{A}$
are calculated by the implicit restart Arnoldi method, which is implemented in ARPACK (Sorensen et al. Reference Sorensen, Lehoucq, Yang and Maschhoff1996). The inversion step used in Arnoldi is achieved by the lower–upper decomposition implemented in SuperLU (Li et al. Reference Li, Demmel, Gilbert, Grigori, Shao and Yamazaki1999). Note that a positive growth rate (
$\omega _i\gt 0$
) indicates a globally unstable flow, while a negative growth rate (
$\omega _i\lt 0$
) signifies a stable flow.
2.3. Resolvent analysis
When the flow is unable to support a global instability (i.e.
$\omega _i \lt 0$
), it can amplify continuous external forcing to realize the transition process. To describe this amplification behaviour, a small-amplitude forcing term is added to (2.5),
\begin{equation} \frac {\partial {\boldsymbol{U}^{\boldsymbol{\prime}}}}{\partial t}=\boldsymbol{A}{\boldsymbol{U}^{\boldsymbol{\prime}}}+\boldsymbol{B}{\boldsymbol{f'}}, \end{equation}
where
$\boldsymbol{B}$
is the matrix that specifies the location and amplitude of the forcing. Assuming that the forcing is periodic both in time and in the spanwise direction, it also can be expressed in the modal form as
\begin{equation} {\boldsymbol{f'}}(x,y,z,t)=\boldsymbol{\hat {f}}(x,y)\exp (\mathrm{i}{\beta }z-\mathrm{i}{\omega }_rt). \end{equation}
After a sufficiently long time, the forced solution exhibits the same wavenumber
$\beta$
and frequency
$\omega _{r}$
as forcing
\begin{equation} {\boldsymbol{U}^{\boldsymbol{\prime}}}(x,y,z,t)=\boldsymbol{\hat {U}}(x,y)\exp (\mathrm{i}{\beta }z-\mathrm{i}{\omega }_rt). \end{equation}
Substituting (2.9) and (2.10) into (2.8) gives
\begin{equation} \hat {\boldsymbol{U}}={\boldsymbol{R}}\boldsymbol{B}\boldsymbol{\hat {f}},\quad {\boldsymbol{R}}=-(\boldsymbol{A}+\mathrm{i}\omega _{r}\boldsymbol{I})^{-1}, \end{equation}
where
$\boldsymbol{I}$
is the identity matrix and
$\boldsymbol{R}$
is the resolvent matrix that determines the relationship between forcing and its linear response. To quantify the maximum amplification, we define the optimal gain as the ratio of the energy-based response norm to the energy-based forcing norm,
\begin{equation} \sigma ^2(\beta ,\omega _r)=\max _{\boldsymbol{f}}\frac {\big\|\boldsymbol{\hat {U}} \big \|_E}{\big \|\boldsymbol{\hat {f}} \big\|_E}. \end{equation}
Chu energy (Chu Reference Chu1965) is utilized to quantify the energy-based norm as
\begin{equation} \big \|\boldsymbol{\hat {U}} \big \|_{E}=\boldsymbol{\hat {U}}^{*}\boldsymbol{M}\boldsymbol{\hat {U}}, \end{equation}
where
$(\boldsymbol{\cdot })^*$
denotes complex-conjugate transpose and
$\boldsymbol{M}$
is the weight matrix, as detailed in Bugeat et al. (Reference Bugeat, Chassaing, Robinet and Sagaut2019). Since matrix
$\boldsymbol{M}$
is invertible, (2.12) is simplified to
\begin{equation} \boldsymbol{B}^*\boldsymbol{M}^{-1}\boldsymbol{R}^*\boldsymbol{MRB\hat {f}}=\sigma ^2\boldsymbol{\hat {f}}. \end{equation}
The optimization problem is then transformed into an eigenvalue problem where the optimal gain is the largest eigenvalue and the optimal forcing
$\boldsymbol{\hat {f}}$
is the corresponding eigenvector. As stated in § 2.2, the eigenvalue is solved using ARPACK (Sorensen et al. Reference Sorensen, Lehoucq, Yang and Maschhoff1996).
2.4. Linear stability analysis
In this paper, LST is used to validate the resolvent analysis. It is a local method and suitable for weakly non-parallel flows. The small perturbation is assumed to be periodic in normal and spanwise directions, as well as in time, expressed in the modal form as
\begin{equation} \boldsymbol{\phi }(x,y,z,t)=\boldsymbol{\hat {\phi }}(y)\exp [\mathrm{i}({{\alpha }_r+\mathrm{i}\alpha }_i)x+\mathrm{i}{\beta }z-\mathrm{i}{\omega }t]. \end{equation}
Here,
$\boldsymbol{\hat {\phi }}(y)=({\hat {u}},{\hat {v}},{\hat {p}},{\hat {T}},{\hat {w}})^T$
represents the eigenfunction of perturbation,
$\beta$
is the spanwise wavenumber,
$\omega$
is the angular frequency. Here
${\alpha }_{r}$
and
${\alpha }_i$
denote the streamwise wavenumber and the spatial growth rate, respectively. Substituting (2.15) into (2.4) yields the following ordinary differential equations:
\begin{equation} \left(\boldsymbol{A_l}\frac {{\rm d}^2}{{\rm d}y^2}+\boldsymbol{B_l}\frac {\rm d}{{\rm d}y}+\boldsymbol{C_l}\right)\boldsymbol{\hat {\phi }}(y)=0. \end{equation}
The matrices
$\boldsymbol{A_l}$
,
$\boldsymbol{B_l}$
, and
$\boldsymbol{C_l}$
are determined by the base flow values, with their specific forms described by Malik (Reference Malik1990). To incorporate the curvature effect, a coordinate transformation is applied using a scale factor defined as
$h_1=1+Ky$
, where
$K$
represents the curvature. This transformation introduces curvature-dependent terms into the differential operators along the wall-normal direction
$y$
. The specific coefficients, including the curvature effect, are comprehensively documented in the Appendix of Ren & Fu (Reference Ren and Fu2014). The boundary conditions are given by
\begin{align}&\,\, \hat {u}=\hat {v}=\hat {T}=\hat {w}=0,y=0; \end{align}
\begin{align}& {\hat {u}}=\hat {v}=\hat {T}=\hat {w}=0,y\rightarrow \infty . \end{align}
The discretization methods for (2.16) include both global and local approaches. The global method uses the single-domain Chebyshev spectral collocation technique, which provides the complete spectrum of eigenvalues, while the local method employs the fourth-order compact difference scheme to refine the eigenvalues obtained from the global method. For detailed information on the LST method, see Guo et al. (Reference Guo, Gao, Jiang and Lee2020) and Guo, Hao & Wen (Reference Guo, Hao and Wen2025).
3. Computational details
3.1. Configurations and flow conditions
The flow conditions are taken from an oblique SWBLI experiment by Willems et al. (Reference Willems, Gülhan and Steelant2015):
$M_{\infty }=6$
,
$p_{\infty }=305$
$\mathrm{Pa}$
,
$T_{\infty }=73$
$\mathrm{K}$
and
$Re_{\infty }=3\times 10^6$
$\mathrm{m^{-1}}$
. The inviscid shock incident position is at
$L=240$
$\mathrm{mm}$
, which is the corner position for the compression corner case in this paper as well. The computational domains and flow structure schematics are depicted in figure 2. When the deflection angle
$\theta$
or ramp angle
$\alpha$
is sufficiently large, a separation bubble is generated due to the adverse pressure gradient. Meanwhile, a separation shock and a reattachment shock are formed as a result of the flow deflection. It is evident that the flow structures of these two configurations are similar, except for the shape of the separation bubble and the expansion waves at the apex of the separation bubble in the incident shock case. The deflection angles
$\theta$
and ramp angles
$\alpha$
are chosen such that they correspond to the same pressure rise. Specifically, the deflection angles for the incident shock flows are selected first, and the corresponding ramp angles are then calculated using the oblique shock relations.
Computational domains and schematics of (a) incident shock flow and (b) compression corner flow.
3.2. Boundary conditions and grids
The incident shock is generated by segmenting the inflow conditions. On the left-hand boundary, below the shock position, the inflow condition aligns with the free stream condition outlined in § 3.1. Above this position and the upper boundary, the flow condition behind the shock is applied, which is calculated using the oblique shock relations. The right-hand boundary is set as an outflow condition achieved by simple extrapolation. The lower boundary is a no-slip and isothermal wall with
$T_w=344$
$\mathrm{K}$
. For the compression ramp, the left-hand and upper boundaries are set to be the free stream condition as normal.
The grids used for GSA consist of
$700\times 300$
cells in the streamwise and normal directions for both incident shock and compression ramp cases. The grid height of the first layer in the normal direction is
$5\times 10^{-6}$
m, resulting in a grid Reynolds number of the order of 1. In resolvent analysis, it is essential to have sufficiently dense streamwise points to accurately capture Mack’s second mode. To satisfy this requirement, there are at least 20 grid points per wavelength of Mack’s second mode. The height of the first normal layer remains
$5\times 10^{-6}$
m. Consequently, the grid for resolvent analysis comprises
$1100\times 250$
cells in the streamwise and normal directions. Grid independence is verified in Appendix A.
4. Results
4.1. Base flows
The characteristics of the 2-D base flows are presented in this section. The skin pressure coefficient
$C_{\!p}$
and skin friction coefficient
$C_{\!f}$
, are defined as
\begin{equation} C_{p}=\frac {p_{w}}{0.5\rho _{\infty }u_{\infty }^{2}},\quad C_f=\frac {\tau _w}{0.5\rho _\infty u_\infty ^2}, \end{equation}
where
$p_w$
and
$\tau _w$
denote the pressure and shear stress on the wall, respectively. Figure 3 depicts the distributions of
$C_{\!p}$
and
$C_{\!f}$
for incident shock and compression corner flows at various angles, where
$0^{\mathrm{\circ }}$
represents the flat plate case for comparison. When a shock incident on the flat plate or generated by the ramp, it induces a pressure rise and imposes an adverse pressure gradient on the boundary layer. In the
$\theta =2^{\mathrm{\circ }}$
incident shock case and
$\alpha =4^{\mathrm{\circ }}$
compression corner case, the adverse pressure gradient is insufficient to produce a separation bubble. When the shock becomes stronger with increasing angles, the adverse pressure gradient becomes large enough to cause separation (
$\theta =4^{\mathrm{\circ }}$
and
$\alpha =8.1^{\mathrm{\circ }}$
). Consequently, separation shock and reattachment shock are induced. Therefore, the single pressure rise is replaced by two rises around separation and reattachment positions. In cases of large separation (
$\theta \gt 4^{\mathrm{\circ }}$
and
$\alpha \gt 8.1^{\mathrm{\circ }}$
), a pressure plateau appears in the middle of the separation region. It should be noted that the discussion regarding the relationship between separation states and angles is specific to the current flow conditions and cannot be directly generalized.
The equivalent ramp angles with the same pressure rise are calculated using oblique shock relations, which neglect the viscous effect. Nevertheless, the numerical pressure rises between the two configurations are close. If the pressure rise is quantified as
\begin{equation} p_r=(p_2-p_1)/p_1, \end{equation}
where
$p_1$
and
$p_2$
represent the pressures before and after the interaction region. The differences of
$p_{r}$
between incident shock and compression corner cases at various angles are all less than
$4\, \%$
. Despite nearly identical pressure rise, the separation length of the incident shock flow is smaller than that of the compression corner flow. A detailed comparison of the separation points, reattachment positions and separation lengths are provided in table 1.
Distributions of (a)
$C_{\!p}$
for incident shock flows, (b)
$C_{\!p}$
for compression corner flows, (c)
$C_{\!f}$
for incident shock flows and (d)
$C_{\!f}$
for compression corner flows. Circle symbols, separation and reattachment locations; horizontal lines, zero skin friction.
Comparison of the separation between compression corner flows and incident shock flows.
The numerical schlieren of incident shock and compression corner flows are presented in figure 4, where the module of density gradient
$|\boldsymbol{\nabla }\rho |$
is defined as
\begin{equation} |\boldsymbol{\nabla }\rho |=\sqrt {\left(\frac {\partial {\rho }}{\partial {x}}\right)^2+\left(\frac {\partial {\rho }}{\partial {y}}\right)^2}. \end{equation}
The flow structure of incident shock flow is slightly more complex compared with that of compression corner flow. In addition to the leading-edge shock, separation shock and reattachment shock, the incident shock flow features expansion waves at the apex of the separation bubble. Moreover, shock–shock interactions occur between the leading-edge shock and incident shock, as well as between the incident shock and the separation shock. Notably, the actual shock incident position is shifted forward relative to the inviscid incident position due to the thickening and separation of the boundary layer. As a result, the separation region of incident shock flow is located upstream compared with the compression corner flow. This distinction is also listed in table 1.
Numerical schlieren of (a) incident shock flow at
$\theta =5^{\mathrm{\circ }}$
and (b) compression corner flow at
$\alpha =10.2^{\mathrm{\circ }}$
.
4.2. Global stability analysis
The separation bubble exhibits an intrinsic instability known as global instability when the strength of SWBLI is sufficiently high. The ramp or deflection angle that determines the onset of global instability is referred to as the critical angle. Beyond the critical angle, boundary layer flow becomes globally unstable. As stated by Song & Hao (Reference Song and Hao2023), for the same pressure rise, incident shock flow exhibits greater stability compared with the equivalent compression corner flow. This discrepancy in stability boundaries can be attributed to several factors from the perspective of base flows. Firstly, the length of separation of incident shock flow is shorter than that of the compression corner flow. Another factor is the presence of expansion waves in the incident shock case, which stabilize the flow. Additionally, the upstream shift of the separation bubble in incident shock flow results in a smaller local Reynolds number, further enhancing stability.
Based on the global stability criteria established by Hao et al. (Reference Hao, Cao, Wen and Olivier2021) and Song & Hao (Reference Song and Hao2023), the critical ramp angle of the compression corner case and critical deflection angle of incident shock case are
$\alpha =10.4^{\mathrm{\circ }}$
and
$\theta =5.8^{\mathrm{\circ }}$
, respectively. To validate the critical angles and explore the characteristics of global instability, GSA is conducted at various angles. The variations of growth rates with spanwise wavenumbers shown in figures 5(a) and 5(b) indicate that the numerical critical ramp angle
$\alpha$
for the compression corner is between
$10.2^{\mathrm{\circ }}$
and
$11^{\mathrm{\circ }}$
, while for the incident shock flow, the numerical critical deflection angle
$\theta$
is between
$5.5^{\mathrm{\circ }}$
and
$6^{\mathrm{\circ }}$
, validating the accuracy of the established criteria. When the angle slightly exceeds the critical angle (compression corner (CC)
$11^{\mathrm{\circ }}$
and incident shock (IS)
$6^{\mathrm{\circ }}$
), only a stationary instability emerges.
As the angle increases further, the separated flow becomes more unstable and contains more unstable modes. For the compression corner at
$\alpha =12^{\mathrm{\circ }}$
, an oscillatory but stable mode appears at large spanwise wavenumbers, referred to as the ‘S-O mode’ in figure 5(a). The dashed line connecting the ‘S mode’ and ‘S-O mode’ indicates that these two modes originate from the same mode essentially. Similarly, an additional oscillatory mode appears in incident shock flow at
$\theta =6.5^{\mathrm{\circ }}$
. This mode is unstable and does not belong to the branch of the primary stationary mode. When the spanwise wavenumber reaches
${\beta }L=75$
, a pair of oscillatory conjugate modes transforms into two stationary modes. The resulting stationary mode with a larger growth rate is shown in figure 5(b), labelled as ‘O-S mode’. The dashed line connecting ‘O-S mode’ and ‘O mode’ suggests that they originate from the same family of instabilities.
Variations of eigenvalues of the most unstable modes with spanwise wavenumbers for (a) compression corner flow and (b) incident shock flow. Here S mode, stationary mode; O mode, oscillatory mode; S-O mode, stationary to oscillatory mode; O-S mode, oscillatory to stationary mode.
The eigenfunctions of the most unstable modes are depicted in figure 6. These eigenfunctions are predominantly concentrated within the separation bubble, with the signs of spanwise perturbations alternating between positive and negative. Between figures 6(d) and 6(e), the eigenfunctions remain similar although the unstable mode has transferred from oscillatory into stationary. The reconstructed 3-D eigenfunctions are presented in Appendix B. Another notable scenario in figure 5 is the increase in growth rate as the spanwise wavenumber approaches zero. The eigenfunctions for the
$\alpha =10.2^{\mathrm{\circ }}$
compression corner flow and
$\theta =5.5^{\mathrm{\circ }}$
incident shock flow at
${\beta }L=0$
are shown in figure 7. These eigenfunctions are primarily located in the shear layer of separation bubble and extend downstream slightly.
Real parts of eigenfunction
$w'$
for the most unstable modes for
$12^{\mathrm{\circ }}$
compression corner flow at (a)
${\beta }L=22$
(S mode) and (b)
${\beta }L=58$
(S-O mode), and for
$6.5^{\mathrm{\circ }}$
incident shock flow at (c)
${\beta }L=34$
(S mode), (d)
${\beta }L=72$
(O mode), and (e)
${\beta }L=79$
(O-S mode). Here black lines are for dividing streamlines.
Real parts of eigenfunction
$u'$
for the least stable modes at
${\beta }L=0$
for (a)
$10.2^{\mathrm{\circ }}$
compression corner flow and (b)
$5.5^{\mathrm{\circ }}$
incident shock flow. Here black lines are for streamlines at the edge of the boundary layer.
4.3. Resolvent analysis
For globally stable flows, resolvent analyses are conducted to investigate how the separation bubble amplifies external forcing. Four deflection angles are selected for incident shock flows:
$0^{\mathrm{\circ }}$
(flat plate),
$2^{\mathrm{\circ }}$
(no separation),
$4^{\mathrm{\circ }}$
(small separation) and
$5^{\mathrm{\circ }}$
(large separation), respectively. Their compression corner counterparts with the same pressure rises are considered as well, with ramp angles determined as
$0^{\mathrm{\circ }}$
,
$4^{\mathrm{\circ }}$
,
$8.1^{\mathrm{\circ }}$
and
$10.2^{\mathrm{\circ }}$
using oblique shock relations. The forcing is fixed at
$x_f/L=0.083$
for both cases, which extends from the wall to the free stream. The results of resolvent analysis are not sensitive to the forcing location, as validated in Appendix C. For simplicity,
$\omega$
is used to denote the angular frequency, rather than
$\omega _{r}$
.
4.3.1. Optimal gain
The optimal gains in a wide space of spanwise wavenumbers and angular frequencies for compression corner cases at various ramp angles are depicted in figure 8. In the flat plate case shown in figure 8(a), the maximum optimal gain is located at
${\beta }L=96$
and
${\omega }L/u_{\infty }=24$
(
$f=16.4$
kHz), which is supposed to be the oblique first mode. Another moderate optimal gain is observed in the low-frequency region at
${\omega }L/u_{\infty }\lt 2.4$
(
$f\lt 1.64$
kHz), spanning a wide range of spanwise wavenumbers of
$50\lt {\beta }L\lt 180$
. It is related to the streamwise streaks in the form of counterrotating streamwise vortices. In the high-frequency region, a local optimal gain appears at
${\omega }L/u_{\infty }=120$
(
$f=81.9$
kHz) and
${\beta }L=0$
, corresponding to the second mode. As the ramp angle increases to
$\alpha =4^{\mathrm{\circ }}$
, which is insufficient to induce separation, the gain contour is shown in figure 8(b). Compared with the flat plate, streaks, the first mode and the second mode are all amplified, with streaks exhibiting the most significant enhancement. The peak optimal gain of streaks shifts to a larger wavenumber around
${\beta }L=200$
, and the range of wavenumbers broadens as well. When the ramp angle further increases to
$\alpha =8.1^{\mathrm{\circ }}$
, a small separation is induced. Figure 8(c) illustrates that streaks continue to experience significant amplification. The first mode is also amplified in the presence of separation. The gain of the second mode remains almost unchanged while the frequency increases. These amplification trends persist in the large separation case at
$\alpha =10.2^{\mathrm{\circ }}$
, as shown in figure 8(d), though the second mode appears weaker compared with the
$\alpha =8.1^{\mathrm{\circ }}$
case.
The optimal gain contours for compression corner flows in the space of spanwise wavenumbers and angular frequencies at (a)
$\alpha =0^{\mathrm{\circ }}$
, (b)
$\alpha =4^{\mathrm{\circ }}$
, (c)
$\alpha =8.1^{\mathrm{\circ }}$
and (d)
$\alpha =10.2^{\mathrm{\circ }}$
.
The optimal gain contours for incident shock flows in the space of spanwise wavenumbers and angular frequencies at (a)
$\theta =0^{\mathrm{\circ }}$
, (b)
$\theta =2^{\mathrm{\circ }}$
, (c)
$\theta =4^{\mathrm{\circ }}$
and (d)
$\theta =5^{\mathrm{\circ }}$
.
The optimal gain contours of incident shock cases are depicted in figure 9. When a shock wave impinges on the flat plate without causing separation (
$\theta =2^{\mathrm{\circ }}$
), the resulting gain contour is shown in figure 9(b). The first mode is slightly amplified, while the gain of streaks is weakened. Notably, the second mode experiences an amplification by an order of magnitude. When the shock wave induces separation, the scenario is similar to that of compression corner flow. In the small separation case at
$\theta =4^{\mathrm{\circ }}$
(see figure 9
c), both streaks and the first mode are enhanced, with the most significant amplification in streaks at larger wavenumber and over a broader range. Conversely, the second mode is diminished slightly. In the large separation case at
$\theta =5^{\mathrm{\circ }}$
(see figure 9
d), the amplification of streaks intensifies further, with the maximum gain shifting to
${\beta }L=312$
. The first mode becomes increasingly difficult to distinguish due to the expanded range of wavenumbers associated with streaks. Meanwhile, the second mode remains slightly diminished, with its frequency increasing as the deflection angle grows.
To validate the dominance of optimal responses, suboptimal gains are calculated for the
$10.2^{\mathrm{\circ }}$
compression corner and
$5^{\mathrm{\circ }}$
incident shock cases. Figure 10 presents the suboptimal gain along with gain separation contours (the difference between the optimal and suboptimal gains). It is evident that, in regions of significant amplification, the suboptimal gains in both cases are at least three orders of magnitude lower than the optimal gains, confirming the predominance of the optimal response. Additionally, for the compression corner flow at
$\alpha =10.2^{\mathrm{\circ }}$
, the first mode is more clearly identifiable in the gain separation contour (figure 10
b), located around
${\beta }L=120$
and
${\omega }L/u_{\infty }=16.8$
.
Contours of (a) suboptimal gain and (b) gain separation between optimal and suboptimal gains for
$\alpha =10.2^{\mathrm{\circ }}$
compression corner flow; (c) suboptimal gain and (b) gain separation between optimal and suboptimal gains for
$\theta =5^{\mathrm{\circ }}$
incident shock flow.
The suboptimal results confirm that the analysis of optimal responses is sufficient for the present cases. By comparing the optimal gains of compression corner and incident shock flows, several common characteristics emerge in the amplification of external forcing by SWBLIs. Firstly, streaks undergo the most significant amplification, with the peak gain shifting to larger wavenumbers and the range of wavenumbers being widened. Secondly, the first mode is enhanced as well, though not as strongly as the streaks, and its frequency and wavenumber remain almost unchanged. Regarding the second mode, it is amplified when the shock wave is insufficient to cause separation while exhibiting insignificant growth in the presence of separation. Additionally, the frequency of the second mode increases slightly with the shock strength. The main difference between the two configurations is that, the amplifications in compression corner flows are generally stronger than those in incident shock flows, which can be attributed to the stabilizing effect of expansion waves in incident shock flows. It is worth noting that in the
$\theta =2^{\mathrm{\circ }}$
incident shock case, streaks are not amplified and the range of wavenumbers shrinks. The reason for this abnormal phenomenon will be explained later.
4.3.2. Amplification of streaks
Figures 11 and 12 present the optimal responses of streaks with the maximum gain for compression corner and incident shock flows at various angles. Figures 11(a–d) and 12(a–d) and figures 11 (e–h) and 12 (e–h) correspond to the spanwise velocity and temperature distributions, respectively. The frequency is fixed at
$f=1.6$
Hz and the results are consistent at lower frequencies. Within the boundary layer, spanwise velocity disturbances of streaks alternate between positive and negative signs in the wall-normal direction, manifesting as streamwise vortices in 3-D flows. In the flat plate boundary layer, the spanwise velocity of streaks undergoes rapid transient growth near the forcing location and gradually decays downstream, whereas the temperature response continues to increase. In compression corner flows, streaks begin to grow at the corner when no separation occurs (figure 11
b,f); whereas in the presence of separation (figure 11
c,d,g,h), streaks grow around separation and reattachment points, accelerating notably after reattachment.
The above patterns also hold in incident shock flows. The crucial difference between the optimal responses of these two configurations lies in the presence of expansion waves. Specifically, in the case without separation (figure 12 b,f), streaks are disrupted by expansion waves, which explains the diminishment of streaks in the weak incident shock case as observed in figure 9(b). This effect will be further observed quantitatively through Chu energy distributions. Additionally, the sign of the spanwise disturbance changes after shock incidence, and the disturbance is radiated outwards by the expansion waves (figures 12 b–d). Similar to compression corner flows, the streaks in incident shock flows are significantly amplified after reattachment.
Optimal responses of streaks with maximum gain for compression corner flows at (a,e)
$\alpha =0^{\mathrm{\circ }}$
, (b,f)
$\alpha =4^{\mathrm{\circ }}$
, (c,g)
$\alpha =8.1^{\mathrm{\circ }}$
and (d,h)
$\alpha =10.2^{\mathrm{\circ }}$
; (a–d) spanwise velocity, (e–h) temperature. Here open circles are for separation and reattachment points.
Optimal responses of streaks with maximum gain for incident shock flows at (a,e)
$\theta =0^{\mathrm{\circ }}$
, (b,f)
$\theta =2^{\mathrm{\circ }}$
, (c,g)
$\theta =4^{\mathrm{\circ }}$
and (d,h)
$\theta =5^{\mathrm{\circ }}$
; (a–d) spanwise velocity; (e–h) temperature. Here open circles are for separation and reattachment points.
To investigate the growth mechanism, Chu energy is integrated along the wall-normal direction. Figure 13 shows the distributions of Chu energy along the
$x$
direction for streaks with maximum gain at various angles, together with the flat plate results for comparison. As stated by Hao et al. (Reference Hao, Cao, Guo and Wen2023), the amplification of streaks in SWBLIs is attributed to the Görtler mechanism. To further explore this, the curvatures of streamlines at the edge of the boundary layers are also plotted in figure 13. The curvature based on the streamline coordinates is calculated using the following equation:
\begin{equation} K=\frac {{\rm d}^2y/{\rm d}x^2}{\left ( 1+\left ( {\rm d}y/{\rm d}x \right ) ^2 \right ) ^{3/2}}. \end{equation}
As mentioned earlier, streaks originate from transient growth near the leading edge. In compression corner flows (figure 13
a–c), streaks are amplified in regions of significant concave curvatures. For example, in the absence of separation, the large curvatures are found at the corner; in the presence of separation, the large curvatures arise from the separation and reattachment points. The amplification of streaks accompanies these regions of large curvatures. In the middle of the separation, the growth of streaks slows down due to the small curvature of the relatively flat shear layer, which is particularly evident in large separation case at
$\alpha =10.2^{\mathrm{\circ }}$
(figure 13
c).
The distributions of Chu energy for the corresponding incident shock flows are shown in figure 13(d–f). Unlike in compression corner flows, the curvature distributions in incident shock flows exhibit three distinct peaks. The first and third peaks are associated with the separation and reattachment, respectively. A large negative curvature occurs at the apex of separation bubble or the shock incident position, representing the convex curvature (cannot be seen in figure 13). The second peak appears just after the apex of the separation bubble, which can be observed in the streamlines shown in figure 6(c–e) and figure 7(b). This peak is presumed to result from the combined effects of the incident shock and the expansion waves, which cause the streamlines to deflect and create a concave curvature region. Notably, this intermediate peak does not appear at Mach 2.15 (Song & Hao Reference Song and Hao2023), suggesting that it is not a universal feature of all incident shock flows. Further investigation is required to fully understand the specific conditions that lead to the formation of this second curvature peak. In incident shock flows, Chu energy is amplified around the separation point but is interrupted by expansion waves, leading to a drop in energy. After the second peak, rather than around the third peak at reattachment, the energy grows dramatically again, driven by the large curvature in this region.
Distributions of integrated Chu energy for the most amplified streaks at various angles with the curvatures of streamlines at the edge of the boundary layers: (a–c), compression corner flow; (d–f), incident shock flow. Here dash–dotted lines, integrated Chu energy of the flat plate case; vertical dashed lines, separation and reattachment locations; blue dashed lines, growth rates of LST.
The blue dashed lines in figure 13 represent the growth rates of the most unstable stationary Görtler mode calculated by LST considering the curvature effect. For compression corner flows, there is excellent agreement between the growth rates obtained from LST and the slopes of Chu energy derived from resolvent analysis. This consistency confirms the convective-type nature of the streaks’ growth at large curvatures. Figure 14 lists the eigenfunctions for
$\alpha =10.2^{\mathrm{\circ }}$
compression corner flow (strong interaction case) at various streamwise locations, corresponding to the locations of three growth rates in figure 13(c). The eigenfunctions obtained from resolvent analysis are extracted at the same streamwise locations. The horizontal and vertical axes represent the normalized disturbance amplitude and wall-normal coordinates, respectively. The wall-normal shapes of the eigenfunction from resolvent analysis and LST are almost identical, especially for the density disturbance. The peaks of disturbances are located near the edge of the shear layer. Discrepancies observed outside the boundary layer in figure 14(a,c) are caused by the influence of separation and reattachment shocks. The discrepancy in velocity near the wall in figure 14(b) can be attributed to non-parallelism and component-type instability.
For incident shock flows (figure 13
d–f), there is generally good agreement in the growth rates, except in regions following the expansion waves. In these regions, LST underestimates the growth rates for separation cases compared with resolvent analysis. Eigenfunctions for a weak interaction case of
$\theta =4^{\mathrm{\circ }}$
incident shock flow are plotted in figure 15. Same with figure 14, the subfigures correspond to different streamwise locations of growth rates shown in figure 13(e). The eigenfunctions in figure 15(a,c) show good alignment. However, the discrepancy appears after expansion waves, shown in figure 15(b). It is observed that significant density disturbances emerge after expansion waves (indicated by the horizontal dashed line) in resolvent analysis results, which are not captured by LST. This observation likely explains why the growth rates obtained from LST are smaller than those in resolvent analysis after expansion waves.
Comparison of eigenfunctions of streaks for
$\alpha =10.2^{\mathrm{\circ }}$
compression corner flow with
${\beta }L=288$
at (a)
$x/L=0.7$
, (b)
$x/L=1.05$
and (c)
$x/L=1.2$
. Here horizontal dashed lines are for positions of separation and reattachment shocks.
Comparison of eigenfunctions of streaks for
$\theta =4^{\mathrm{\circ }}$
incident shock flow with
${\beta }L=288$
at (a)
$x/L=0.8$
, (b)
$x/L=1$
and (c)
$x/L=1.2$
. Here horizontal dashed lines are for positions of separation and reattachment shocks in (a) and (c); position of expansion waves in (b).
Therefore, for streaks, the amplification is attributed to the Görtler mechanism induced by large concave curvatures, aligning with the conclusion drawn for Mach 7.7 compression corner flow by Hao et al. (Reference Hao, Cao, Guo and Wen2023). The LST further validates the results of resolvent analysis but tends to underestimate the growth rates for incident shock flows near expansion waves. The overall amplification of streaks in incident shock flows is less than that in compression corner flows, primarily due to the interruption caused by the expansion waves.
4.3.3. Amplification of the first mode
The optimal responses of the first mode with maximum gain are illustrated in figure 16. The first mode is primarily concentrated in the shear layer. Similar to streaks, the first mode undergoes amplification around the separation and reattachment points, and experiences attenuation at the apex of the separation bubble in incident shock flow due to the influence of expansion waves.
Optimal responses of the first mode with maximum gain for compression corner incident shock flows at (a)
$\alpha =0^{\mathrm{\circ }}$
, (b)
$\alpha =4^{\mathrm{\circ }}$
, (c)
$\alpha =8.1^{\mathrm{\circ }}$
and (d)
$\alpha =10.2^{\mathrm{\circ }}$
; and for incident shock flow at (e)
$\theta =0^{\mathrm{\circ }}$
, (f)
$\theta =2^{\mathrm{\circ }}$
, (g)
$\theta =4^{\mathrm{\circ }}$
and (h)
$\theta =5^{\mathrm{\circ }}$
. Here open circles are for separation and reattachment points.
Distributions of integrated Chu energy for the most amplified first mode at various angles with the curvatures of streamlines at the edge of the boundary layers: (a–c), compression corner flow; (d–f), incident shock flow. Here dash–dotted lines, Chu energy of the flat plate case; vertical dashed lines, separation and reattachment locations; blue dashed lines, growth rates of LST.
Comparison of eigenfunctions of the first mode for
$\alpha =10.2^{\mathrm{\circ }}$
compression corner flow with
${\beta }L=120$
and
${\omega }L/u_{\infty }=16.8$
at (a)
$x/L=0.7$
, (b)
$x/L=1.05$
and (c)
$x/L=1.2$
. Here horizontal dashed lines are for positions of separation and reattachment shocks.
Comparison of eigenfunctions of the first mode for
$\theta =4^{\mathrm{\circ }}$
incident shock flow with
${\beta }L=96$
and
${\omega }L/u_{\infty }=24$
at (a)
$x/L=0.8$
, (b)
$x/L=1$
and (c)
$x/L=1.2$
. Here horizontal dashed lines are for positions of separation and reattachment shocks in (a) and (c); position of expansion waves in (b).
Figure 17 plots the distributions of integrated Chu energy for the first mode with the most significant amplification. In the flat plate case (depicted by dash–dotted lines), the first mode undergoes transient growth, while in regions with large curvatures induced by SWBLI, it exhibits almost linear growth. In flows with separation (figure 17 b,c), the linear growth slows down in the middle of the separating bubble in compression corner flow; whereas in incident shock flow, the growth is interrupted by expansion waves and restores near reattachment points. The amplification pattern of the first mode closely resembles that of streaks, indicating a potential connection to the Görtler mechanism as well. This hypothesis will be further investigated in the subsequent section.
Additionally, LST is performed to facilitate a comparison with resolvent analysis results, with the growth rates shown by blue dashed lines in figure 17. The eigenfunction comparisons for
$\alpha =10.2^{\mathrm{\circ }}$
compression corner flow and
$\theta =4^{\mathrm{\circ }}$
incident shock flow are provided in figures 18 and 19, respectively. For compression corner flows, the similarity of growth rates and eigenfunctions between resolvent analysis and LST confirms that the amplification of the first mode in the separation region is also driven by convective-type instability. For incident shock flows, the growth rates from LST agree well with that of resolvent analysis when away from expansion waves. In contrast, discrepancies appear near expansion waves shown in figure 17(e,f). Reflected in the eigenfunctions, when there is no interference from expansion waves, the eigenfunctions obtained by LST and resolvent analysis are highly consistent, as shown in figure 19(a). With the effects of expansion waves (figure 19
b,c), the density disturbances from the resolvent analysis exhibit oscillations, which are not captured by LST. This observation is consistent with the effect of expansion waves on streaks discussed in § 4.3.2.
Optimal responses of the second mode with maximum gain for compression corner flows at (a)
$\alpha =0^{\mathrm{\circ }}$
, (b)
$\alpha =4^{\mathrm{\circ }}$
, (c)
$\alpha =8.1^{\mathrm{\circ }}$
and (d)
$\alpha =10.2^{\mathrm{\circ }}$
; and for incident shock flows at (e)
$\theta =0^{\mathrm{\circ }}$
, (f)
$\theta =2^{\mathrm{\circ }}$
, (g)
$\theta =4^{\mathrm{\circ }}$
and (h)
$\theta =5^{\mathrm{\circ }}$
. Here open circles are for separation and reattachment points.
Therefore, the amplification of the first mode, similar to the streaks, is driven by large curvatures. The agreement between LST and resolvent analysis results confirms that this amplification originates from convective-type instability. However, in incident shock flows, expansion waves suppress the growth of the first mode and impact the accuracy of LST predictions for the growth rate.
4.3.4. Amplification of the second mode
The optimal responses of the second mode with the maximum gain at various angles are presented in figure 20. In the flat plate case, the second mode grows gradually and exhibits the typical ‘two-cell’ structure within the boundary layer. In compression corner flow without separation, the second mode does not experience dramatic amplification at the corner and the disturbances are radiated along the shock. Interestingly, in flows with separation, higher Mack’s modes appear within the separation bubble, as indicated by additional peaks in pressure disturbances. This phenomenon has also been observed by Balakumar et al. (Reference Balakumar, Zhao and Atkins2005) using LST. The optimal responses in resolvent analysis provide a more intuitive visualization of these higher Mack’s modes. Specifically, in figure 20(c), the third mode appears and alternates with the second mode. In the large separation case shown in figure 20(d), even the fourth and fifth modes appear. These different Mack’s modes are separated by regions of attenuation, with only the second mode persisting downstream of the separation. The characteristics of Mack’s modes for incident shock flows in figure 20(e–f) are similar, with more evident radiation caused by the shock and expansion waves.
The integrated Chu energy shown in figure 21 reveals that within the separation bubble, the second mode cannot sustain continuous growth and instead exhibits oscillatory development. In some cases, the boundary layer downstream of reattachment supports the continuous growth of the second mode. This downstream growth can be attributed to a preferred boundary layer thickness of the second mode, which is approximately half of its streamwise wavelength. Therefore, the second mode thrives in a boundary layer with a matching thickness downstream of reattachment while demonstrating minimal growth within the separation bubble. Higher Mack’s modes, which appear within the separation bubble, exhibit unsustainable growth and do not persist downstream.
Distributions of integrated Chu energy for the most amplified Mack’s modes at various angles with the curvatures of streamlines at the edge of the boundary layers: (a–c), compression corner flow; (d–f), incident shock flow. Here dash–dotted lines, Chu energy of the flat plate case; vertical dashed lines, separation and reattachment locations; blue dashed lines, growth rates of LST.
A comparison between resolvent analysis and LST is also conducted, with the slopes of blue dashed lines in figure 21 indicating the growth rates from LST, and eigenfunctions shown in figures 22 and 23. Outside the separation, the growth rates are accurately captured by LST. The pressure disturbances exhibit a typical second mode pattern, characterized by a first peak at the wall and a second peak near the boundary layer edge (figures 22 a and 23 a). Within the separation bubble, the growth rates obtained by LST differ from the slopes of Chu energy in resolvent analysis. Nevertheless, the eigenfunctions near the wall from LST closely resemble those from resolvent analysis. Additionally, the third and fourth peaks of LST eigenfunctions within the separation (figures 22 b and 23 b) validate the presence of higher Mack’s modes observed in optimal responses from resolvent analysis.
Comparison of eigenfunctions of the second mode for
$\alpha =10.2^{\mathrm{\circ }}$
compression corner flow with
${\omega }L/u_{\infty }=168$
at (a)
$x/L=0.55$
, (b)
$x/L=1.2$
and (c)
$x/L=1.45$
. Here horizontal dashed lines are for positions of reattachment shocks.
Comparison of eigenfunctions of the second mode for
$\theta =4^{\mathrm{\circ }}$
incident shock flow with
${\omega }L/u_{\infty }=168$
at (a)
$x/L=0.6$
, (b)
$x/L=0.85$
and (c)
$x/L=1.2$
. Here horizontal dashed lines are for positions of separation and reattachment shocks.
4.3.5. Rounded compression corner
To validate whether the amplification of streaks and the first mode is motivated by the Görtler mechanism, rather than the inherent properties of separation bubble, the sharp compression corner with a ramp angle of
${\alpha }=10.2^{\mathrm{\circ }}$
is modified by rounding the corner to prevent separation, with a radius of
$r/L=4.2$
. A comparison of skin pressure coefficient
$C_{\!p}$
and skin friction coefficient
$C_{\!f}$
between rounded and sharp corner cases is presented in figure 24. The results indicate that the rounded corner case achieves an equivalent pressure rise to that of the sharp corner case while effectively eliminating the separation.
Distributions of (a) skin pressure coefficient and (b) skin friction coefficient for sharp corner and rounded corner cases. Here open circles, separation and reattachment points of sharp corner case; solid circles, tangency points of rounded corner case.
The gain contours for rounded compression corner in the space of spanwise wavenumber and angular frequency of (a) optimal gain, (b) suboptimal gain and (c) gain separation.
Resolvent analysis is performed on the rounded corner case, with the resulting optimal gain contour depicted in figure 25(a). The optimal streaks occur within a low-frequency region around a spanwise wavenumber of
${\beta }L=288$
, which is larger than
${\beta }L=240$
observed in the sharp corner case (refer to figure 8
d). A clear peak at
${\omega }L/u_{\infty }=144$
signifies the presence of the second mode, which exhibits a stronger amplification compared with the sharp corner case. As in the sharp corner scenario, the identification of the first mode is challenging due to the dominant amplification of streaks. Figure 25(b) illustrates the suboptimal gain contour, which is significantly weaker than the optimal gain and can thus be neglected. Nevertheless, the gain separation between optimal and suboptimal gains, as shown in figure 25(c), effectively highlights the first mode by eliminating the influence of dominant streaks. Specifically, the first mode is located at
${\beta }L=144$
and
${\omega }L/u_{\infty }=12$
.
Optimal responses with maximum gain for rounded compression corner case: (a) streaks; (b) first mode; (c) second mode. Here open circles are for tangency points.
Figure 26 presents the optimal responses for rounded compression corner flow, with the corresponding development of Chu energy illustrated in figure 27. In the rounded corner case, both streaks and the first mode exhibit nearly exponential growth along the circular arc. The elimination of the separation bubble leads to the disappearance of higher Mack’s modes. Chu energy for the remaining second mode undergoes continuous growth along the circular arc due to the constant boundary layer thickness, and then shows oscillation downstream.
Distributions of integrated Chu energy for the most amplified (a) streaks, (b) the first mode and (c) the second mode for rounded compression corner flow with the curvature of streamline at the edge of the boundary layers. Here dash–dotted lines, Chu energy of the flat plate case; vertical dashed lines, tangency points of rounded corner case; blue dashed lines, growth rates of LST.
4.3.6. Discussion of the curvature-amplified first mode
The similar amplification trends observed between streaks and the first mode in the rounded corner case further verify the hypothesis proposed in § 4.3.3, that the amplification of the first mode in SWBLIs may also be driven by the Görtler mechanism. To delve deeper into the underlying motivation, a detailed analysis is conducted for the first mode of rounded corner case using LST. The angular frequency and spanwise wavenumber are consistent with that of the first mode exhibiting optimal gain captured by resolvent analysis. The position where LST is performed is fixed at
$x/L=1$
, unless otherwise specified.
Eigenvalue trajectory at
${\omega }L/u_{\infty }=12$
and
${\beta }L=144$
with decreasing curvature for rounded corner flow at
$x/L=1$
.
Firstly, the curvature is reduced artificially with a scale ratio, which is incorporated into the LST framework by modifying the parameter
$K$
as mentioned in § 2.4. The eigenvalue trajectory with decreasing curvatures is shown in figure 28. The growth rate
$\alpha _i$
of the discrete unstable mode decreases continuously with the curvature
$K$
. When the curvature is close to zero,
$\alpha _i$
is less than zero, indicating that the mode is stable. The continuous decrease in growth rate demonstrates that this unstable mode is curvature-amplified.
Furthermore, to explore the origin of this curvature-amplified first mode, LST without the curvature effect is carried out. As shown in figure 29(a), when the angular frequency decreases from
${\omega }L/u_{\infty }=12$
to
${\omega }L/u_{\infty }=0$
, this discrete mode gradually approaches the entropy–vorticity spectrum. This observation complements the receptivity studies of Ma & Zhong (Reference Ma and Zhong2003) and Fedorov (Reference Fedorov2011), which suggested that the first mode in supersonic flows originates from the acoustic spectrum. This discrepancy is presumed to be due to the difference between oblique waves and planar waves. Hence, the eigenvalue trajectory is examined with decreasing spanwise wavenumbers, as depicted in figure 29(b). Initially, as the spanwise wavenumber decreases, the growth rate of the discrete mode increases. This trend arises because the most unstable oblique mode at
${\omega }L/u_{\infty }=12$
shifts to a smaller wavenumber as the curvature effect diminishes in the disturbance equations. When the spanwise wavenumber reaches
${\beta }L=0$
, the planar waves are observed to originate from the slow acoustic spectrum, consistent with Ma & Zhong (Reference Ma and Zhong2003). The above results suggest a potential distinction between the oblique and planar first modes: while the planar first mode originates from the acoustic spectrum, the oblique first mode may arise from the entropy–vorticity spectrum, consistent with the understanding that the first mode is the compressible counterpart of Tollmien–Schlichting mode.
Eigenvalue trajectory of the first mode at (a)
${\beta }L=144$
and
$0\leqslant {\omega }L/u_{\infty }\leqslant 12$
and (b)
${\omega }L/u_{\infty }=12$
and
$0\leqslant {\beta }L\leqslant 144$
for rounded corner flow at
$x/L=1$
.
Additionally, as observed in figure 29(a), a discrete unstable mode still persists at zero frequency. This occurs because the curvature effect remains embedded in the base flow profile, despite being reduced to zero in the disturbance equations of LST. It is supposed that this discrete mode can fully emerge into the entropy–vorticity spectrum only when the curvature effect in the base flow is also eliminated. Therefore, the receptivity process is repeated at
$x/L=0.55$
, a location on upstream flat plate. As shown in the eigenvalue trajectory in figure 30(a), the discrete mode completely progressively merges into the entropy–vorticity spectrum as the frequency decreases. In fact, when
${\omega }L/u_{\infty }\lt 2.4$
, the discrete mode becomes nearly indistinguishable from the entropy–vorticity spectrum. The eigenvalue trajectory with respect to spanwise wavenumber at
$x/L=0.55$
is plotted in figure 30(b). Same with the case at
$x/L=1$
, the discrete mode originates from the slow acoustic spectrum as
${\beta }L$
approaches zero. The LST results for
$x/L=0.55$
further validate the classification between oblique and planar waves. Based on this clarification, referring to these oblique waves as the ‘first mode’ may not be entirely precise because they originate from the entropy–vorticity spectrum. In the subsequent discussion, we use the term ‘oblique mode’ to describe what was previously referred to as the ‘first mode’.
Eigenvalue trajectory of the first mode at (a)
${\beta }L=144$
and
$0\leqslant {\omega }L/u_{\infty }\leqslant 12$
and (b)
${\omega }L/u_{\infty }=12$
and
$0\leqslant {\beta }L\leqslant 144$
for rounded corner flow at
$x/L=0.55$
.
Some representative conditions are selected to compare the eigenfunctions, shown in figures 31 and 32. Specifically, at
$x/L=1$
, the following five conditions are selected: (i)
${\beta }L=144$
and
${\omega }L/u_{\infty }=12$
with
$K$
; (ii)
${\beta }L=144$
and
${\omega }L/u_{\infty }=12$
with
$0K$
; (iii)
${\beta }L=72$
and
${\omega }L/u_{\infty }=12$
with
$0K$
; (iv)
${\beta }L=144$
and
${\omega }L/u_{\infty }=0$
with
$K$
; (v)
${\beta }L=0$
and
${\omega }L/u_{\infty }=12$
with
$0K$
. The eigenfunctions at
${\beta }L=144$
and
${\omega }L/u_{\infty }=12$
remain extremely similar regardless of the presence of curvature. As the spanwise wavenumber decreases to
${\beta }L=72$
, where the largest growth rate appears in figure 29(b), the eigenfunction maintains its shape. Furthermore, the eigenfunction associated with the stationary Görtler instability (
$K$
,
${\beta }L=144$
,
${\omega }L/u_{\infty }=0$
) still exhibits a strong resemblance to the first three conditions. However, when
${\beta }L=0$
, the eigenfunction develops a pronounced free stream signature, likely due to its proximity to the slow acoustic spectrum. This noticeable transformation of eigenfunctions with decreasing spanwise wavenumber may indicate a shift in the underlying instability mechanism. For
$x/L=0.55$
, the representative conditions for eigenfunction comparison are as follows: (i)
${\beta }L=144$
and
${\omega }L/u_{\infty }=12$
with
$0K$
; (ii)
${\beta }L=48$
and
${\omega }L/u_{\infty }=12$
with
$0K$
(corresponding to the largest growth rate in figure 30
b); (iii)
${\beta }L=0$
and
${\omega }L/u_{\infty }=12$
with
$0K$
. Again, The eigenfunctions of the oblique modes exhibit strong similarity, whereas differ significantly from that of the planar mode.
Eigenfunctions of typical conditions at
$x/L=0.55$
: (a) density and (b) streamwise velocity.
Eigenfunctions of typical conditions at
$x/L=0.55$
: (a) density and (b) streamwise velocity.
The variation of phase speed with spanwise wavenumber.
To further examine the hypothesis regarding the distinction between oblique and planar waves, figure 33 presents the variations of phase speed with spanwise wavenumber at
$x/L=1$
and
$x/L=0.55$
, along with the phase speed of slow acoustic waves obtained from the eigenvalue spectrum of LST. For small spanwise wavenumbers, the phase speed of the oblique mode closely follows the trend of slow acoustic waves. As the spanwise wavenumber increases, the phase speed reaches a minimum value in the range of
$0.7{\sim }0.8 u_\infty$
. With a further increase in wavenumber, the phase speed gradually rises and approaches
$u_\infty$
, aligning with the phase speed of the entropy–vorticity spectrum. This trend may indicate a transfer in instability mechanisms, potentially shifting from acoustic waves to Görtler instability.
Therefore, the oblique mode observed in the present study appears to be essentially equivalent to unsteady Görtler instability, both of which originate from the entropy–vorticity waves. This explains why the oblique mode undergoes significant amplification in the presence of concave curvatures, exhibiting nearly linear growth as streaks.
It should be noted that this classification is derived from the present case, which requires further investigations across more cases to validate this conclusion. Moreover, in 3-D flows, the flow modulation in spanwise direction can alter the amplification mechanism of the first mode. Therefore, in the future work, 2-D LST and 3-D PSE will be conducted on a 3-D baseflow distorted by streaks under the weakly nonlinear framework, and how the first mode and streaks interact with each other will be examined through nonlinear resolvent analysis.
5. Conclusions
This paper provides a comprehensive investigation of instabilities in SWBLIs at Mach 6, with a main focus on the amplification mechanisms revealed by resolvent analysis and LST. Additionally, it explores the characteristics of global instability. Both compression corner and shock incidence are examined to identify the commonalities and differences between these two configurations.
First, the global stability boundary is determined using GSA. For compression corner flow, this boundary lies between ramp angles of
$10.2^{\mathrm{\circ }}$
and
$11^{\mathrm{\circ }}$
, while for incident shock flow, it is between deflection angles of
$5.5^{\mathrm{\circ }}$
and
$6^{\mathrm{\circ }}$
. The ranges confirm the accuracy of the critical criteria proposed by Hao et al. (Reference Hao, Cao, Wen and Olivier2021) and Song & Hao (Reference Song and Hao2023). Beyond these boundaries, a stationary unstable mode appears. In the
$12^{\mathrm{\circ }}$
compression corner flow, the stationary unstable mode transitions into an oscillatory mode but remains stable at large wavenumbers. In the
$6.5^{\mathrm{\circ }}$
incident shock flow, an additional oscillatory mode emerges, which switches to stationary at large wavenumbers.
Under the global stability boundary, resolvent analysis is conducted to explore the instabilities associated with upstream perturbation amplification. Four distinct flow conditions are considered: flat plate, no separation, small separation and large separation. Equivalent cases between the incident shock and the compression corner flows share the same pressure rise. The resolvent analysis results demonstrate that streaks are the most amplified instability, followed by the oblique mode. While the second mode exhibits insignificant growth and higher Mack’s modes are observed in the presence of the separation bubble. The main difference in incident shock flows compared with compression corner flows is that expansion waves suppress the growth of the instability, resulting in a slightly lower optimal gain for the incident shock flow.
The skin friction coefficient distributions of different grids for (a) GSA and (b) resolvent analysis.
The eigenvalue spectra from GSA of different grids for (a)
$6^{\mathrm{\circ }}$
incident shock flow at
${\beta }L=33.6$
and (b)
$12^{\mathrm{\circ }}$
compression ramp flow at
${\beta }L=21.6$
.
The optimal gains from resolvent analysis of different grids for (a)
$5^{\mathrm{\circ }}$
incident shock flow and (b)
$10.2^{\mathrm{\circ }}$
compression ramp flow at different frequencies. Here symbols are for fine grid.
Isosurfaces of the real parts of eigenfunctions
$u'$
for
$10.2^{\mathrm{\circ }}$
compression ramp flow at (a)
$\beta L=22$
(S mode) and (b)
$\beta L=58$
(S-O mode), and for
$6.5^{\mathrm{\circ }}$
incident shock flow at (c)
$\beta L=34$
(S mode), (d)
$\beta L=72$
(O mode) and (e)
$\beta L=79$
(O-S mode).
The amplification of streaks is attributed to the Görtler instability, which exhibits linear growth at large concave curvatures. The growth pattern of the oblique mode is similar to that of streaks, indicating that the oblique mode is also driven by the Görtler instability. The second mode shows no obvious growth within the separation bubble, while experiences amplification after reattachment if the boundary layer thickness matches. In the rounded compression corner case, it is validated that the growth of both streaks and the oblique mode is independent of separation. Receptivity analysis and phase speed variation using LST further suggest that the oblique mode may originate from the entropy–vorticity mode, reinforcing the conclusion that the oblique mode’s amplification is driven by the Görtler instability.
Therefore, this paper has revealed the complete instability characteristics of SWBLIs at Mach 6. Future work will focus on the transition process in SWBLIs, with attention to how various instabilities and their interactions contribute to the onset of transition.
Acknowledgements
The authors gratefully acknowledge Dr P. Guo (Department of Aeronautical and Aviation Engineering, The Hong Kong Polytechnic University) for generously providing the LST code used in this study.
Funding
This work is supported by the National Natural Science Foundation of China (no. 12102377) and the Hong Kong Research Grants Council (nos. 15204322 and 25203721).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Grid independent verification
The gird cells used in this paper are detailed in § 3.2, with
$700\times 300$
cells for GSA and
$1100\times 250$
for resolvent analysis. For gird independent verification, fine grids consisting of
$1000\times 400$
cells are used for GSA and
$1300\times 350$
cells for resolvent analysis. The
$6^{\mathrm{\circ }}$
incident shock and
$12^{\mathrm{\circ }}$
compression corner cases are examined for GSA, while grid independence for resolvent analysis is verified by the
$5^{\mathrm{\circ }}$
incident shock and
$10.2^{\mathrm{\circ }}$
compression corner cases. The skin friction coefficients of the base flows are shown in figure 34, demonstrating that the grids have achieved convergence for the base flows. The eigenvalue spectra of different grids at the wavenumber with the maximum growth rate are presented in figure 35. The growth rates are very close, particularly for the most unstable modes. Resolvent analysis is performed on fine grids at the frequencies and wavenumbers corresponding to streaks, the first mode and the second mode. The optimal gains compared in figure 36 indicate that the grid resolution is sufficiently accurate, with only minor deviations for the high-frequency, which are within reasonable limits.
The optimal gains as a function of spanwise wavenumber with different forcing locations for
$\alpha =8.1^{\mathrm{\circ }}$
compression corner case.
Appendix B. The 3-D eigenfunctions of GSA
Figure 37 presents the reconstructed 3-D eigenfunctions of GSA, corresponding to the 2-D eigenfunctions shown in figure 6.
Appendix C. Forcing location validation
The comparison of optimal gains with different forcing locations for
$\alpha =8.1^{\mathrm{\circ }}$
compression corner case is shown in figure 38. In this test case, the forcing is applied slightly upstream at
$x/L = 0.05$
. The selected frequencies correspond to those of the streaks, the first mode and the second mode. It can be observed that the optimal gains are insensitive to the forcing location.
Open access
