3.1. Governing equations

The compressible Navier–Stokes equations in operator form are given as

(3.1) \begin{equation} \frac {\partial \boldsymbol{U}}{\partial t}=\boldsymbol{\mathcal{N}}(\boldsymbol{U}), \end{equation}

where $ \boldsymbol{U}$ is the vector of conservative variables, and $ \boldsymbol{\mathcal{N}}$ denotes the nonlinear Navier–Stokes operator. For the RANS equations, $ \boldsymbol{U}$ is the vector of Favre-averaged variables, and $ \boldsymbol{\mathcal{N}}$ is the nonlinear RANS operator. In the case of LES, the governing equations are the Favre-filtered Navier–Stokes equations, with $ \boldsymbol{U}$ representing the vector of Favre-filtered variables, and $ \boldsymbol{\mathcal{N}}$ representing the nonlinear LES operator. The Favre filter is defined as $ \langle f \rangle =\overline {\rho f} /\bar {\rho }$ , where an overline denotes the spatial filter (Garnier, Adams & Sagaut Reference Garnier, Adams and Sagaut2009).

The working fluid is assumed to be a perfect gas with Prandtl number $ Pr=0.72$ and specific heat ratio 1.4. The molecular viscosity $ \mu$ is calculated from Sutherland’s law, and the turbulent Prandtl number $ Pr_{t}$ is chosen to be 0.9.

3.2. The RANS solver

The Boussinesq assumption is employed in the RANS equations to model the Reynolds stresses, with the eddy viscosity $ \mu _{t}$ obtained from the Spalart–Allmaras turbulence model (Spalart & Allmaras Reference Spalart and Allmaras1992). An in-house finite-volume solver called PHAROS (Hao, Wang & Lee Reference Hao, Wang and Lee2016; Hao & Wen Reference Hao and Wen2020) is used to obtain RANS base flows. The modified Steger–Warming scheme (MacCormack Reference MacCormack2014) is applied in computing the inviscid fluxes, and the viscous fluxes are solved using a second-order central difference scheme. The implicit line relaxation method proposed by Wright, Candler & Bose (Reference Wright, Candler and Bose1998) is employed for pseudo-time marching. The computation domain and boundary conditions are also marked in figure 1. The far-field boundaries follow the freestream conditions. The supersonic outflow boundary uses a simple extrapolation, and an isothermal no-slip condition with fixed temperature $ 275.4$ K is applied to the wall.

3.3. The LES solver

Implicit LES are also conducted using PHAROS. The low-dissipative sixth-order kinetic preserving scheme (Pirozzoli Reference Pirozzoli2010) combining the AUSM $ ^{+}$ -up scheme (Liou Reference Liou2006) with the fifth-order weighted essentially non-oscillatory reconstruction (Jiang & Shu Reference Jiang and Shu1996) is used to solve the inviscid fluxes, switched by the Jameson sensor (Jameson, Schmidt & Turkel Reference Jameson, Schmidt and Turkel1981). Time integration is performed using a third-order low-storage Runge–Kutta scheme, with time step $ 9$ ns. The total simulation time is approximately 10 ms. Flow samples are collected for duration 8.1 ms ( $ 2200 \delta /u_{\infty }$ ).

The computational domain for LES is reduced appropriately, as depicted in figure 1. The spanwise width of the domain is set as $ 6\delta$ . The inflow boundary employs the extended digital filter technology (Touber & Sandham Reference Touber and Sandham2009; Ceci et al. Reference Ceci, Palumbo, Larsson and Pirozzoli2022) and a one-dimensional profile extracted from the RANS base flow, as marked at station $ x_{in}$ at $ (-22.5\delta ,0)$ in figure 1. The profile is chosen to ensure that the boundary-layer thickness $ \delta$ at station $ x_{0}$ is 2.27 mm after a turbulence recovery process. Supersonic outflow boundary conditions are specified at the upper and right boundaries. The grid is progressively coarsened at these boundaries, combined with a sponge zone of thickness approximately $ 3\delta$ near the upper boundary to damp the flow variables to the freestream values to avoid any reflections (Mani Reference Mani2012). Periodic boundary conditions are used in the spanwise direction. The wall is assumed to be no-slip and isothermal with $ T_{w}=275.4$ K. The grid measures $1423 \times181\times301$ in the streamwise, wall-normal and spanwise directions, respectively. The mesh is uniform in the streamwise and spanwise directions, with grid resolutions $ \Delta x^{+}\approx 19$ , $\Delta {y}_{w}^{+}\approx 1.0$ and $ \Delta z^{+}\approx 9$ at the reference station $ x_{0}$ .

3.4. Resolvent analysis

Resolvent analysis is based on the scale separation assumption, which posits a sufficient scale gap between large-scale coherent structures and turbulent fluctuations. Globally stable flows can be modelled as a noise amplifier (Huerre & Monkewitz Reference Huerre and Monkewitz1990) and investigated using resolvent analysis (Sipp & Marquet Reference Sipp and Marquet2013). To study the linear forced dynamics of a noise amplifier flow, a small-amplitude forcing $ \boldsymbol{f}^{\prime}$ is introduced into the linearised Navier–Stokes equations, given as

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

where $ \boldsymbol{\mathcal{A}}$ is the linearised operator evaluated using a base flow from RANS results or LES Favre-averaged results, and $ \boldsymbol{U}^{\prime}$ is the perturbation term around the base flow. The operator $ \boldsymbol{\mathcal{B}}$ is to constrain the forcing field to a localised region (Bugeat et al. Reference Bugeat, Chassaing, Robinet and Sagaut2019). Both $ \boldsymbol{U}^{\prime}$ and $ \boldsymbol{f}^{\prime}$ are assumed to be harmonic in time and the spanwise direction, given as

(3.3) \begin{align} \boldsymbol{U}^{\prime}(x,y,z,t) & =\hat {\boldsymbol{U}}(x,y)\exp(i\beta z-i{\omega }_{r}t), \nonumber\\ \boldsymbol{f}^{\prime}(x,y,z,t) & =\,\hat {\boldsymbol{\!\!f}}(x,y)\exp(i\beta z-i{\omega }_{r}t), \end{align}

where $ \beta$ is the spanwise wavenumber, and $ \omega _{r}$ is the angular frequency. Substituting (3.3) into (3.2) gives

(3.4) \begin{equation} -i{\omega }_{r}\hat {\boldsymbol{U}}=\boldsymbol{\mathcal{A}\hat {U}}+\boldsymbol{\mathcal{B}\hat {f}}. \end{equation}

Introducing the resolvent operator $ \boldsymbol{\mathcal{R}}={(-{i\omega }_{r}\boldsymbol{\mathcal{I}}-\boldsymbol{\mathcal{A}})}^{-1}$ and the identity operator $ \boldsymbol{\mathcal{I}}$ , (3.4) can be written as a relationship between the constrained forcing and its linear response, given by

(3.5) \begin{equation} \hat{\boldsymbol{U}}=\boldsymbol{\mathcal{R}}\boldsymbol{\mathcal{B}}\,\,\hat {\boldsymbol{\!\!f}}. \end{equation}

To seek the most energetic amplification between the forcing and the response, the optimal gain $ \sigma$ is defined as

(3.6) \begin{equation} {\sigma }^{2}({\omega }_{r},\beta )=\max _{\hat {f}}{\frac {{\| \boldsymbol{\hat {U}}\| }_{E}^{2} }{{\| \boldsymbol{{\hat {f}}}\| }_{E}^{2}}}. \end{equation}

The Chu energy norm (Chu Reference Chu1965) is applied to the response and the forcing to evaluate their respective energies. As discussed by Bugeat et al. (Reference Bugeat, Chassaing, Robinet and Sagaut2019), (3.6) can be converted to an eigenvalue problem. The problem is discretised similarly to the flow solver, except that a central scheme is used to solve the inviscid fluxes in smooth regions detected by a shock sensor (Hendrickson, Kartha & Candler Reference Hendrickson, Kartha and Candler2018). At the left-hand and upper boundaries, all perturbations are set to zero. At the right-hand boundary, a simple extrapolation is used. The wall conditions are $ u^{\prime}=v^{\prime}=w^{\prime}=0$ , $ T^{\prime}=0$ and $ \partial p^{\prime}/\partial n =0$ . Sponge zones with thickness $ 5\delta$ are implemented to the left-hand, right-hand and upper boundaries to gradually damp the perturbations to zero (Mani Reference Mani2012). Finally, the implicitly restarted Arnoldi method implemented in ARPACK (Lehoucq, Sorensen & Yang 1996), interfaced with the MUMPS (Amestoy et al. Reference Amestoy, Duff, L’Excellent and Koster2001) open library, is used to solve the eigenvalue problems.

It should be noted that the resolvent analysis solver adapted is for laminar flows (Hao et al. Reference Hao, Cao, Guo and Wen2023), while incorporating the effective viscosity $ \mu _{eff}=\mu +\mu _{t}$ and effective heat conductivity $ k_{eff}=\mu /Pr+\mu _{t}/Pr_{t}$ to model turbulent fluctuations (Fan et al. Reference Fan, Kozul, Li and Sandberg2024). The explicit expressions of the operators $ \boldsymbol{\mathcal{A}}$ , $ \boldsymbol{\mathcal{B}}$ and $ \boldsymbol{\mathcal{R}}$ in laminar flows can be found in Bugeat et al. (Reference Bugeat, Chassaing, Robinet and Sagaut2019). In this study, the perturbation of $ \mu _{t}$ is not considered, meaning that $ \mu _{t}$ is assumed to be frozen everywhere, consistent with Cura et al. (Reference Cura, Hanifi, Cavalieri and Weiss2024).

3.5. Local stability analysis

In local stability analysis (LSA), the perturbation $ \boldsymbol{U}^{\prime}$ is expressed as

(3.7) \begin{equation} \boldsymbol{U}^{\prime}(x,y,z,t)=\hat {\boldsymbol{U}}_{1\text{-}D}(y)\exp\left [i(\alpha _{r}+i\alpha _{i})x+i\beta z-i{\omega }_{r}t \right ], \end{equation}

where $ \hat {\boldsymbol{U}}_{1\text{-}D}(y)$ is the one-dimensional eigenfunction, $ \alpha _{r}$ is the real streamwise wavenumber, and $ \alpha _{i}$ is the spatial growth rate. A negative value of $ \alpha _{i}$ indicates that the corresponding mode is unstable. Notably, the effective viscosity $ {\mu }_{eff}$ and the concave curvature $ K$ (Ren & Fu Reference Ren and Fu2014) are adopted to conduct LSA. Substituting (3.7) into the linearised Navier–Stokes equations yields an eigenvalue problem, which is solved using the Eigen (Guennebaud, Jacob & Reference Guennebaud2010) library at given values of $ \beta$ and $ {\omega }_{r}$ . The boundary conditions are consistent with those in § 3.2.

3.6. Spectral proper orthogonal decomposition

Spectral proper orthogonal decomposition (Lumley Reference Lumley1967; Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018) extracts coherent structures from the perspective of energy from flow data. It seeks a set of orthogonal bases to represent the overall power spectra in time and space optimally. Hence SPOD is ideal for identifying coherent structures in turbulent flows.

The key point in SPOD is to accurately estimate the cross-spectral density (CSD) tensor from a time series (Towne et al. Reference Towne, Schmidt and Colonius2018). The CSD matrix at each frequency is estimated by

(3.8) \begin{equation} \boldsymbol{\hat {S}}(\omega _{k})=\boldsymbol{\hat {Q}}(\omega _{k})\boldsymbol{\cdot} {\boldsymbol{\hat {Q}}^{*}}(\omega _{k}), \end{equation}

where $ \boldsymbol{\hat {Q}}(\omega _{k})$ is the data matrix at $ k$ th frequency $ \omega _{k}$ , and $ (\cdot )^{*}$ is the Hermitian transpose. The problem is transformed to a series of eigenvalue problems at frequency $ \omega _{k}$ , given as

(3.9) \begin{equation} \boldsymbol{\hat {S}W\hat {\Psi }}=\boldsymbol{\hat {\Psi }\hat {\Lambda }}, \end{equation}

or via the method of snapshots,

(3.10) \begin{equation} \boldsymbol{\hat {Q}}^{*} \boldsymbol{W} \hat{\boldsymbol{\!Q}} \hat {\boldsymbol{\Theta}}=\boldsymbol{\hat {\Theta }} \hat {\boldsymbol{\Lambda} },\quad\boldsymbol{\hat {\Psi }} =\boldsymbol{\hat {Q}}\hat {\boldsymbol{\Theta }}, \end{equation}

where $ \boldsymbol{W}$ is the weight matrix, $ \boldsymbol{\hat {\!\Lambda } }$ gives the eigenvalues, and the SPOD modes $ \boldsymbol{\hat {\Psi }}$ are recovered using the eigenvectors $ \boldsymbol{\hat {\Theta }}$ . The Chu energy norm (Chu Reference Chu1965) is commonly used in the weight matrix $ \boldsymbol{W}$ for compressible flows.

Of particular interest is the most energetic SPOD mode, which corresponds to the largest eigenvalue, often referred to as the optimal or leading SPOD mode. A low-rank feature is observed when the optimal eigenvalue is significantly larger than the others. In this scenario, the physical mechanism associated with the optimal mode becomes dominant and prevalent (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018).

4.1. Response to upstream disturbances

Figure 2 presents the Mach number contour of the time- and spanwise-averaged flow obtained from LES. The verification of our LES solver is discussed in Appendix A. In the concave surface region, a series of compression waves is generated due to an adverse pressure gradient. These waves converge downstream of the concave surface, forming a shock wave. Figure 3 shows the distributions of the skin friction coefficient $ {C}_{f}$ and the wall pressure coefficient $ {C}_{p}$ . These coefficients are defined as

(4.1) \begin{equation} {C}_{f}=\frac {2{\tau }_{w}}{{\rho }_{\infty }{u}_{\infty }^{2}},\quad {C}_{p}=\frac {2{p}_{w}}{{\rho }_{\infty }{u}_{\infty }^{2}}, \end{equation}

where $ {\tau }_{w}$ and $ {p}_{w}$ are the wall shear stress and pressure, respectively. The behaviour of $ {C}_{f}$ shows a sudden decrease at the onset of the concave surface, followed by an approximately linear increase across the concave region. Similarly, the pressure coefficient $ {C}_{p}$ displays a corresponding increasing trend throughout this region.

Figure 2.

Contour of Mach number distribution at $ M_{\infty }=2.95$ , $ Re_{\delta }=63\,500$ . The open circle indicates the end of the concave surface.

Figure 3.

Distributions of ( $ a$ ) the skin friction coefficient $ C_{f}$ , and ( $ b$ ) the wall pressure coefficient $ C_{p}$ , at $ r=50\delta$ and $ \varphi =20^\circ$ . Open circles indicate the end of the concave region.

Figure 4.

( $ a$ ) Optimal gains as a function of the spanwise wavenumber $ \beta \delta$ at frequencies $ St=3.6\times {10}^{-4}, 3.6\times {10}^{-3}, 3.6\times {10}^{-2}$ . ( $ b$ ) Optimal gains at $ \beta \delta=2.6$ over a range of frequencies. The black dashed line in ( $ a$ ) means local maxima, while in ( $ b$ ) it denotes the approximate cut-off frequency $ St =0.036$ . The red dashed lines indicate the cut-off frequency.

Resolvent analysis is then performed to study responses of the linear dynamic system (3.5) to upstream forcings, utilising the LES-computed mean flow as the base flow (Touber & Sandham Reference Touber and Sandham2009). The eddy viscosity $ {\mu }_{t}$ in the LES base flow is computed from

(4.2) \begin{equation} {\mu }_{t}={C}_{\mu }\rho \frac {{k}^{2}}{\epsilon }, \end{equation}

where $ {C}_{\mu }=0.09$ , $ \rho$ is the averaged density, and $ \epsilon$ is the rate of dissipation of turbulent kinetic energy (Pope Reference Pope2001). The turbulent kinetic energy $ k$ is defined as $ k=({1}/{2})\times (\langle u^{\prime\prime} u^{\prime\prime}\rangle + \langle v^{\prime\prime} v^{\prime\prime}\rangle + \langle w^{\prime\prime} w^{\prime\prime}\rangle )$ , where $ \langle u^{\prime\prime} u^{\prime\prime}\rangle$ , $ \langle v^{\prime\prime} v^{\prime\prime}\rangle$ and $ \langle w^{\prime\prime} w^{\prime\prime}\rangle$ are the Reynolds stress components. The forcing station is set to $ x=-10\delta$ . It should be noted that the optimal response is insensitive to the forcing location (see Appendix B).

Figure 4( $ a$ ) compares optimal gains at three frequencies $ St=3.6\times {10}^{-4}, 3.6\times {10}^{-3}, 3.6\times {10}^{-2}$ over a range of spanwise wavenumbers $ \beta \delta$ . Here, the frequency is defined as $ St=f\delta /{u}_{\infty }$ . The overlap of optimal gains at $ St=3.6\times {10}^{-4}$ and $ St=3.6\times {10}^{-3}$ indicates that the optimal gains do not change significantly with further reductions in frequency. For $ St=0.036$ , the optimal gains are slightly lower than those at the other two frequencies. Interestingly, all three optimal gain curves peak at the same spanwise wavenumber $ \beta \delta = 2.6$ . Figure 4( $ b$ ) presents the optimal gains at the peak wavenumber $ \beta \delta = 2.6$ as a function of $ St$ . The distribution exhibits characteristics consistent with a first-order low-pass filter (Hunter Reference Hunter2001), demonstrating peak gains in the low-frequency regime $ St\lt 0.01$ . Beyond $ St=0.01$ , a gradual attenuation of optimal gain is observed, transitioning to a linear decay pattern when exceeding $ St=0.04$ . An approximate cut-off frequency $ St=0.036$ is identified through the intersection of the $ St\lt 0.01$ and $ St\gt 0.04$ trends (marked in red) (Hunter Reference Hunter2001). This critical frequency denotes the threshold beyond which the linear system experiences a marked reduction in optimal gains. The selection of this specific cut-off frequency is further substantiated by SPOD analysis results, as detailed in § 4.2.2. The highest frequency considered in this study is $ St=0.072$ , which is at least one order lower than the characteristic frequency of TBLs, thereby satisfying the scale separation assumption.

Figure 5 shows the optimal forcing and corresponding response at $ St=0.036$ with spanwise wavenumber $ \beta \delta=2.6$ . The forcing corresponds to streamwise vortices, since the amplitudes of the wall-normal component $ | v^{\prime} |$ and the spanwise component $ | w^{\prime} |$ are larger than that of the streamwise component $ | u^{\prime} |$ . The lift-up mechanism works to transform the streamwise vortices into streamwise streaks in the flow response (Abreu et al. Reference Abreu, Cavalieri, Schlatter, Vinuesa and Henningson2020; Hao et al. Reference Hao, Cao, Guo and Wen2023). The optimal response of $ u^{\prime}$ displays streamwise streaks within and downstream of the concave region, while $ w^{\prime}$ is characterised by streamwise counter-rotating vortices, as shown in figures 5(b,c). Upstream of the concave region, no apparent coherent structures can be observed. It is evident that concave curvature plays a vital role in forming and strengthening these coherent structures. In the concave region, the curved streamlines excite centrifugal effects. The secondary effect, known as the Görtler instability, becomes the dominant amplification mechanism in forming the streamwise vortices due to the large concave curvature 0.02 (Hoffmann et al. Reference Hoffmann, Muck and Bradshaw1985). Therefore, the streamwise counter-rotating vortices generated in the concave region are identified as Görtler vortices. In addition, these structures are similar in appearance, differing only in their streamwise lengths across the frequency band $St= 0.01{-}0.04$ , which strongly indicates that they are driven by the same physical mechanism.

Figure 5.

( $ a$ ) Amplitude of the optimal forcing. (b,c) Real parts of $ u^{\prime}$ and $ w^{\prime}$ of the optimal response at $ St= 0.036$ and $ \beta \delta = 2.6$ . Open circles indictae the end of the concave region.

Furthermore, one can notice that the streamwise vortices are progressively strengthened in the concave region due to continuous centrifugal effects, resulting in a more pronounced $ u^{\prime}$ and $ w^{\prime}$ in the latter half of the concave surface, as shown in figures 5(b,c). The three-dimensional perturbation field of $ u^{\prime}$ reconstructed using the optimal response at $ St=0.036$ and $ \beta \delta =2.6$ is presented in figure 6. An intuitive comparison of the Görtler structures can be observed at different stations in the figure. At the beginning of the concave surface (station $ \varphi =0^{\circ }$ ), weak spanwise structures are visible. Downstream of this station, the Görtler structures begin to grow. By the end of the concave surface, both the size and strength of the Görtler structures are significantly enhanced compared to values upstream.

Figure 6.

Three-dimensional reconstructed perturbation field of $ u^{\prime}$ using the optimal response at $ St=0.036$ and $ \beta \delta =2.6$ . The stations, from left to right, are $ \varphi =0^{\circ }$ , $ \varphi =5^{\circ }$ , $ \varphi =10^{\circ }$ , $ \varphi =15^{\circ }$ and $ \varphi =20^{\circ }$ , respectively.

Figure 7.

( $ a$ ) Distributions of the Chu energy density integrated in the wall-normal direction at $ St=0.036$ and $ \beta \delta =2.6$ , and the curvature of the streamline passing through $ (x_{0} , 0.3 \delta)$ . ( $ b$ ) The most unstable spatial growth rate at $ St=0.036$ from LSA as a function of spanwise wavenumbers $ \beta \delta$ at stations $ \varphi = 5^{\circ }$ , $ \varphi = 10^{\circ }$ and $ \varphi = 15^{\circ }$ . The black dashed line in ( $ a$ ) indicates the most unstable spatial growth rate predicted by LSA at station $ \varphi = 10^{\circ }$ , while in ( $ b$ ) it means the local maximum $ \beta \delta =2.98$ at $ \varphi = 10^{\circ }$ . The open circle indicates the end of the concave region.

To study the spatial energy growth rate of the optimal response in the concave region, the integrated wall-normal Chu energy density along the streamwise direction at $ St=0.036$ and $ \beta \delta=2.6$ is presented in figure 7( $ a$ ). The Chu energy density increases exponentially in the concave region, and peaks downstream of the region. The concave curvature $ K$ of the streamline passing through $ (x_{0} , 0.3 \delta)$ stays constant in the concave region, indicating uniform centrifugal effects at different streamwise stations. Consequently, the energy of the coherent structures is exponentially amplified in the concave region under consistently strong centrifugal effects. The slope of the black dashed line is predicted by LSA at station $ \varphi = 10^{\circ }$ at $ St=0.036$ and $ \beta \delta= 2.6$ . There is good agreement between the growth rate from LSA and resolvent analysis at the same station, which indicates that the energy growth is mainly attributed to the convective-type non-normality, i.e. the Görtler instability. Figure 7( $ b$ ) depicts the most unstable spatial growth rates at $ St=0.036$ predicted by LSA at three stations as functions of spanwise wavenumbers $ \beta \delta$ . Although the peak wavenumbers of the three lines differ slightly, they are centred around $ \beta \delta= 2.98$ . There is a small discrepancy between the peak wavenumbers predicted by LSA and resolvent analysis. This discrepancy is reasonable, as the wavenumber predicted by LSA is locally optimal, while the global optimal wavenumber from resolvent analysis may vary due to non-parallelism and non-modal effects (Hanifi, Schmid & Henningson Reference Hanifi, Schmid and Henningson1996; Tumin & Reshotko Reference Tumin and Reshotko2003). Figure 8 compares the wall-normal distributions of the amplitudes of $ u^{\prime}$ and $ w^{\prime}$ obtained from two methods at the three stations $ \varphi = 5^{\circ }$ , $ \varphi = 10^{\circ }$ and $ \varphi = 15^{\circ }$ . Here, $ u^{\prime}$ peaks near the wall then decreases linearly, while $ w^{\prime}$ exhibits two peaks: one near the wall, and another near the edge of the TBL. The trends of $ u^{\prime}$ and $ w^{\prime}$ are almost the same at three streamwise stations of the two results, with minor differences in the location of the second peak of $ w^{\prime}$ and the wall-normal range of the large amplitude of $ u^{\prime}$ . It is confirmed that the most amplified disturbances in resolvent analysis are the Görtler vortices.

Figure 8.

Comparisons of wall-normal distributions of streamwise velocity and spanwise velocity perturbations from LSA and resolvent analysis at stations ( $ a$ ) $ \varphi = 5^{\circ }$ , ( $ b$ ) $ \varphi = 10^{\circ }$ and ( $ c$ ) $ \varphi = 15^{\circ }$ , at $ St=0.036$ and $ \beta \delta= 2.6$ . Distributions are normalised by their respective maximum streamwise velocity perturbation amplitudes.

Figure 9.

Iso-surface of streamwise velocity $ u/{u}_{\infty }=0.55$ , coloured by wall-normal distance.

4.2. Validation of predicted coherent structures

4.2.1. Amplification of turbulence intensity

Figure 9 illustrates the instantaneous iso-surface of streamwise velocity $ u/{u}_{\infty }=0.55$ , coloured by wall-normal distance. On the flat plate, only a few near-wall vortices are visible, whereas large-scale vortices are prominent in the concave region. The centrifugal effects in this region continue to drive low-momentum flow outwards, thereby enhancing the momentum exchange between low- and high-momentum flows (Wang et al. Reference Wang, Wang, Sun, Yang, Zhao and Hu2019). Additionally, it is noteworthy that the spanwise scale of the vortices in the concave region gradually increases.

Figure 10.

Wall-normal distributions of ( $ a$ ) turbulent kinetic energy $ k$ , ( $ b$ ) root mean square values RMS $ ^{+}$ , normalised by $ {u}_{\infty }^{2}$ . The black dashed line in ( $ a$ ) denotes $ y{}_{n}/\delta =0.26$ . In (b), the dashed and solid lines correspond to the reference station ( $x_0$ ) and the end of the concave region ( $\varphi=20^{\circ}$ ), respectively.

Figure 10( $ a$ ) shows the wall-normal distributions of the turbulent kinetic energy $ k$ at the reference station $ x_{0}$ and the four different turning angle stations in the concave surface. At station $ x_{0}$ , $ k$ peaks at a near-wall location, primarily due to high- and low-speed streaks. In the concave region, the inner part of $ k$ is slightly enhanced, while the outer part experiences significant amplification. Continuous centrifugal effects amplify turbulent fluctuations, leading to a noticeable increase in the outer part of $ k$ across various stations. A second peak of $ k$ , associated with the large-scale vortices shown in figure 9, occurs at approximately $ y{}_{n}/\delta =0.26$ at the end of the concave surface.

Figure 10( $ b$ ) compares wall-normal distributions of Reynolds stress components $ \langle u^{\prime\prime} u^{\prime\prime}\rangle$ , $ \langle v^{\prime\prime} v^{\prime\prime} \rangle$ , $ \langle w^{\prime\prime} w^{\prime\prime}\rangle$ and $ \langle u^{\prime\prime} v^{\prime\prime}\rangle$ at the reference station $ x_{0}$ and the end of the concave surface $ \varphi =20^{\circ }$ . Compared to the undisturbed TBL, all Reynolds stress components show significant changes at $ \varphi =20^{\circ }$ . The trends of $ \langle v^{\prime\prime} v^{\prime\prime}\rangle$ and $ \langle w^{\prime\prime} w^{\prime\prime}\rangle$ are the same as those at station $ x_{0}$ , although their amplitudes are increased. In contrast, $ \langle u^{\prime\prime} u^{\prime\prime}\rangle$ exhibits a notable discrepancy in the outer layer. After a decrease from the inner peak, $ \langle u^{\prime\prime} u^{\prime\prime}\rangle$ increases again to form a secondary peak influencing $ y{}_{n}/\delta =0.1{-}0.6$ . The outer peak of $ \langle u^{\prime\prime} u^{\prime\prime}\rangle$ is associated with the large-scale structures highlighted in figure 9. Moreover, the peak of $ k$ in the outer layer at the station $ \varphi =20^{\circ }$ is mainly contributed by the streamwise component $ \langle u^{\prime\prime} u^{\prime\prime}\rangle$ . Notably, the near-wall part of the component $ \langle u^{\prime\prime} v^{\prime\prime}\rangle$ changes sign at station $ \varphi =20^{\circ }$ , which is attributed to mathematical contaminations between longitudinal and wall-normal velocities in Cartesian coordinates (Wu et al. Reference Wu, Liang and Zhao2019). Similar distributions in $ k$ and Reynolds stress components were also discussed by Wu et al. (Reference Wu, Liang and Zhao2019) and Sun et al. (Reference Sun, Sandham and Hu2019).

4.2.2. Modal analysis

The POD analysis is performed on $ y$ – $z$ slices in the concave region to study spatial structures of large-scale vortices. Here, the entire state vector $ {q}^{\prime }= [ \rho ^{\prime},u^{\prime},v^{\prime},w^{\prime},T^{\prime} ]$ is used, combined with the Chu energy norm (Chu Reference Chu1965). The fluctuation variables are normalised with the freestream parameters. The convergence analysis of POD modes is discussed in Appendix C.

Figure 11.

Leading POD modes of ( $ a$ ) $ u^{\prime}$ , ( $ b$ ) $ w^{\prime}$ at station $ \varphi =10^{\circ }$ , and ( $ d$ ) $ u^{\prime}$ , ( $ e$ ) $ w^{\prime}$ at station $ \varphi =20^{\circ }$ . (c,f) The mode energy distributions at the two stations.

Figure 11 demonstrates the components $ u^{\prime}$ and $ w^{\prime}$ of the leading POD modes at stations $ \varphi =10^{\circ }$ and $ \varphi =20^{\circ }$ , as well as mode energy distributions. Three pairs of spanwise structures can be observed in figures 11(a,b) at station $ \varphi =10^{\circ }$ . The component $ u^{\prime}$ is characterised by streamwise streaks, and $ w^{\prime}$ features streamwise counter-rotating vortices. The most energetic spatial structures closely resemble the coherent structures predicted by resolvent analysis in figure 6. At station $ \varphi =20^{\circ }$ , the spatial structures behave analogously to those at $ \varphi =10^{\circ }$ . The spanwise scale at both stations is approximately $ {\lambda }_{z}=2\delta$ , leading to three pairs of spanwise structures. Regarding the energy distribution, the energy ratio of the leading POD mode at both stations is minimal, falling below $ 5\, \%$ . Most of the energy is captured by numerous high-order modes. In addition, the sum of eigenvalues at station $ \varphi =20^{\circ }$ is 0.57, compared to 0.24 at station $ \varphi =10^{\circ }$ . The rise of total energy agrees well with the turbulence amplification caused by the centrifugal effects discussed in § 4.2.1.

Nevertheless, it should be emphasised that conventional POD modes may conflate dynamically distinct phenomena when different flow features exhibit comparable energy content at disparate frequencies (Towne et al. Reference Towne, Schmidt and Colonius2018; Mendez, Balabane & Buchlin Reference Mendez, Balabane and Buchlin2019). This multi-modal coupling characteristic fundamentally limits the temporal discrimination capability of energy-based decomposition methods. To overcome this limitation, we use the frequency-resolved method named SPOD to investigate the spatial–temporal features of the coherent structures induced by the concave surface.

Figure 12.

( $ a$ ) Normalised SPOD eigenvalues as functions of $ St$ . The 99 $ \, \%$ confidence bounds for the leading mode are indicated by the green dash-dotted lines. The red dashed line marks the local maximum $ St=0.036$ . ( $ b$ ) Normalised SPOD eigenvalues at $ St=0.036$ .

The SPOD is first performed on the spanwise mid-plane using a temporal segment $ 249.1 \delta /{u}_{\infty }$ . Welch’s method (Welch Reference Welch1967) with 50 $ \, \%$ overlap combining a Hamming window is employed to minimise spectral leakage. Notably, the specific choice of window function shows negligible impact on the results, as equivalent frequency spectra are obtained when using a Hann window. A convergence analysis of the SPOD modes is provided in Appendix C. The spectrum is presented in figure 12( $ a$ ), where the 99 $ \, \%$ bounds indicate the variance of the leading mode. A characteristic frequency $ St=0.036$ of the leading mode is identified. Figure 12( $ b$ ) compares the normalised eigenvalues at $ St=0.036$ . A low-rank feature is observed, with the optimal SPOD mode containing 65 $ \, \%$ of the total flow energy – significantly more energetic than the other suboptimal modes at this frequency.

Figure 13.

Real parts of ( $ a$ ) $ u^{\prime}$ and ( $ b$ ) $ w^{\prime}$ of the optimal SPOD mode at $ St=0.036$ . Open circles indicate the end of the concave surface.

Figures 13(a,b) display the real parts of $ u^{\prime}$ and $ w^{\prime}$ of the optimal SPOD mode at $ St=0.036$ . In response to upstream turbulent forcings, streamwise streaks of $ u^{\prime}$ and $ w^{\prime}$ are observed within and downstream of the concave region. Notably, both $ u^{\prime}$ and $ w^{\prime}$ exhibit stronger intensities in the latter half of the concave surface, indicating an enhancement of the coherent structures in this region. This enhancement process and the coherent structures resemble those obtained from resolvent analysis at the same frequency, as shown in figure 5. To quantitatively compare the leading SPOD mode with the optimal resolvent modes, we calculate the projection coefficient $ \gamma$ using the equation

(4.3) \begin{equation} \gamma =\frac {\left | \| \hat{\boldsymbol{\Psi}}_{1\textit{SPOD}}, \hat{\boldsymbol{U}}_{1\textit{RES}}\|_{E} \right |}{ \| \hat{\boldsymbol{\Psi }}_{1\textit{SPOD}}\|_{E} \times \|\hat{\boldsymbol{U}}_{1\textit{RES}}\|_{E} }, \end{equation}

where $ \hat{\boldsymbol{\Psi}}_{1\textit{SPOD}}$ represents the leading SPOD mode, while $ \hat{\boldsymbol{U}}_{1\textit{RES}}$ denotes the optimal resolvent modes over several spanwise wavenumbers. The notation $ {\| \cdot \|}_{E}$ denotes the Chu norm (Chu Reference Chu1965). A value $ \gamma = 0$ indicates orthogonality between the two modes, while $ \gamma = 1$ signifies perfect alignment (Abreu et al. Reference Abreu, Cavalieri, Schlatter, Vinuesa and Henningson2020; Fan et al. Reference Fan, Kozul, Li and Sandberg2024). Notably, the leading SPOD mode $ \hat{\boldsymbol{\Psi }}_{1\textit{SPOD}}(x,y)$ combines Fourier modes at various wavelengths implicitly, while $ \hat{\boldsymbol{U}}_{1\textit{RES}}(x,y)$ represents structures at a single specified wavelength explicitly. Hence this projection quantifies the structural similarity between the coherent structures captured by the leading SPOD mode and those predicted by resolvent analysis, while simultaneously identifying the dominant wavelength in the SPOD mode through the peak $ \gamma (\beta \delta )$ . Cura et al. (Reference Cura, Hanifi, Cavalieri and Weiss2024) also successfully made the projection between the leading SPOD mode obtained using particle image velocimetry data in a planar and resolvent modes to investigate the low-frequency dynamics of turbulent separation bubbles. Figure 14 shows the projection coefficients $ \gamma$ between the optimal SPOD mode and the predicted resolvent modes at $ St = 0.036$ over various spanwise wavenumbers. The peak alignment, $ \gamma = 0.94$ , is very close to 1, indicating an almost perfect match between the two modes. The optimal response is strongly associated with observed streamwise streaks in SPOD. Moreover, the coefficient reaches its maximum at spanwise wavenumber $ \beta \delta = 2.95$ , which corresponds well with the optimal wavenumber predicted by resolvent analysis and the most unstable wavenumber identified by LSA.

Hence there is a strong correlation between the optimal SPOD mode and the optimal response. The relation between SPOD and resolvent analysis has been discussed by Towne et al. (Reference Towne, Schmidt and Colonius2018), Schmidt et al. (Reference Schmidt, Towne, Rigas, Colonius and Brès2018) and Abreu et al. (Reference Abreu, Cavalieri, Schlatter, Vinuesa and Henningson2020). When forcings are unit-variance white noise, SPOD is equivalent to resolvent analysis. In real turbulent flows, biased and nonlinear forcings (Morra et al. Reference Morra, Nogueira, Cavalieri and Henningson2021) result in differences between SPOD and resolvent modes. However, a close match between the leading SPOD mode and the optimal resolvent mode can be achieved when the resolvent operator identifies a dominant amplification mechanism (Abreu et al. Reference Abreu, Cavalieri, Schlatter, Vinuesa and Henningson2020). In such cases, the leading SPOD modes and the optimal resolvent modes are insensitive to the correlation of nonlinear forcings $ E(\,\,\hat {\boldsymbol{\!\!f}}\,\hat {\boldsymbol{\!\!f}}^{*})$ , where $ E(\cdot )$ is the expectation operator. Therefore, the similarity of mode shapes and the high alignment between the leading SPOD mode and the optimal response suggest that the same physical mechanism, namely the centrifugal effects, dominates the most energetic coherent structures. In other words, the Görtler instability primarily amplifies upstream forcings and forms the Görtler vortices.

Figure 14.

Projection coefficients $ {\gamma }$ between the optimal SPOD mode and resolvent modes over varying spanwise wavenumber at $ St=0.036$ . The black dashed line indicates local maximum $ \beta \delta =2.95$ .

Furthermore, the flow statistics characterised at the peak frequency $ St=0.036$ in the streamwise direction make it reasonable to choose the frequency as the approximate cut-off frequency. The optimal SPOD modes within the frequency band $St= 0.01{-}0.04$ also resemble those obtained from resolvent analysis, with the SPOD mode being more energetic at $ St=0.036$ due to coloured forcings. Consequently, resolvent analysis successfully captures flow structures in the streamwise direction, and serves as an effective tool for studying Görtler instability. The next step is to compare the spanwise scales from both methods, and investigate the formation process of Görtler vortices in the spanwise direction under the influence of centrifugal effects.

The SPOD is then conducted on $ y$ – $ z$ slices in the concave region using the sine-taper SPOD (Yeung & Schmidt Reference Yeung and Schmidt2024). The peak frequencies of the first SPOD modes at different stations do not consistently align with $ St=0.036$ , as shown in table 1. This variation is expected because Görtler vortices are active over a frequency band $St= 0.01{-}0.04$ , as discussed in § 4.1. The coherent structures and spanwise scales of the leading SPOD modes within this frequency band are quite similar. Here, we focus on the coherent structures at $ St= 0.036$ to quantitatively assess their spanwise scales and structures. Figure 15 depicts the real parts of $ u^{\prime}$ and $ w^{\prime}$ of the leading SPOD modes at three stations $ \varphi = {10}^{\circ }$ , $\varphi ={15}^{\circ }$ and $\varphi ={20}^{\circ }$ , at $ St=0.036$ . At station $ \varphi = {5}^{\circ }$ , the Görtler vortices are not fully developed (not shown here). The Görtler vortices emerge at stations $\varphi ={10}^{\circ }$ , where $ u^{\prime}$ displays pairs of streamwise streaks, and $ w^{\prime}$ shows pairs of streamwise counter-rotating vortices. Downstream of this station, the coherent structures at station $\varphi ={15}^{\circ }$ and $\varphi ={20}^{\circ }$ remain unchanged.

Table 1.

The peak frequencies at different stations of the leading SPOD modes.

Figure 15.

Real parts of the optimal SPOD modes of $ u^{\prime}$ and $ w^{\prime}$ at stations (a,d) $ \varphi = {10}^{\circ }$ , (b,e) $ \varphi = {15}^{\circ }$ and (c,f) $ \varphi = {20}^{\circ }$ , at $ St=0.036$ .

The spanwise scales of the Görtler vortices at these stations are approximately consistent, averaging approximately $ 2\delta$ . However, it is noteworthy that some vortices exceed $ 2\delta$ , while others are smaller, suggesting that the spanwise structures exhibit a degree of irregularity. This phenomenon is expected and related to nonlinear effects in turbulent flows. The predicted optimal spanwise wavelength $ 2.4\delta$ by resolvent analysis is slightly larger than the $ 2\delta$ captured by the optimal SPOD mode. The minor discrepancy is attributed to upstream biased forcings. In resolvent analysis, upstream forcings are treated as unbiased, while real turbulent fluctuations are biased and coloured (Morra et al. Reference Morra, Nogueira, Cavalieri and Henningson2021). These coloured forcings tend to concentrate most of the energy in near-wall streaks and very-large-scale motions, explaining the expected $ 2\delta$ patterns observed in the optimal SPOD mode. Overall, resolvent analysis effectively models the Görtler vortices in terms of spatial–temporal features.

Unsteady motions of the Görtler vortices can be recovered using the first five SPOD modes within the frequency band $ St=0.01{-}0.04$ , which contain at least 50 $ \, \%$ of total energy at a given frequency (Nekkanti & Schmidt Reference Nekkanti and Schmidt2021; Yeung & Schmidt Reference Yeung and Schmidt2024). The choice of this frequency band is based on two factors: first, the cut-off frequency band discussed in § 4.1, and second, the dynamic modes strongly related to Görtler vortices at frequencies 0.0115, 0.0172 and 0.0223 found by Tong et al. (Reference Tong, Li, Duan and Yu2017). The unsteady spanwise motions of the Görtler vortices can be observed from time series (see the supplementary movie). The large-scale Görtler vortices undergo a series of unsteady behaviours, including lateral wandering, merging into larger vortices, and splitting into several smaller ones. These merging and splitting behaviours are related to the spanwise motions of vortices, with varying lateral velocities leading to interactions between vortex pairs. No preferred vortex locations are found from the reconstructed data, consistent with experimental findings (Barlow & Johnston Reference Barlow and Johnston1988a ,Reference Barlow and Johnston b ; Schuelein & Trofimov Reference Schuelein and Trofimov2011).

5.1. The RANS-based resolvent analysis

Figure 16.

(a) Normalised optimal gains based on LES and RANS base flows at $ St=0.036$ . ( $ b$ ) Distributions of $ \mu _{t}$ at the reference station $ x_{0}$ from LES and RANS base flows. (c,d) Real parts of $ u^{\prime}$ and $ w^{\prime}$ based on the RANS base flow at $ St=0.036$ and $ \beta \delta =2.75$ . Open circles indicate the ends of the concave surface.

Figure 16( $ a$ ) compares the normalised optimal gains from LES- and RANS-based resolvent analysis at $ St=0.036$ . The RANS results exhibit a marginally higher peak wavenumber $ \beta \delta =2.75$ compared to the LES-based value $ \beta \delta =2.6$ . This minor shift aligns with expectations from eddy viscosity $ {\mu }_{t}$ discrepancies between the methodologies, as discussed in a flat plate case by Cossu, Pujals & Depardon (Reference Cossu, Pujals and Depardon2009). Figure 16( $ b$ ) further reveals that both LES and RANS profiles share similar $ {\mu }_{t}$ trends, peaking near $ y_{n}/\delta =0.5$ , though the RANS-based $ {\mu }_{t}$ magnitude exceeds the LES value by approximately 50 $ \, \%$ . Despite these differences, the optimal responses of $ u^{\prime}$ and $ w^{\prime}$ at $ \beta \delta =2.75$ shown in figures 16(c,d) resemble the optimal responses in figures 5(b,c) and the leading SPOD mode in figures 13(a,b) at the same frequency $ St=0.036$ . Quantitative validation is further supported by a peak projection coefficient $ \gamma =0.9$ between the leading SPOD mode and the RANS-based optimal response at $ \beta \delta =2.95$ . These results collectively indicate that the most amplified dynamics revealed by resolvent analysis is insensitive to the relatively small differences in the $ {\mu }_{t}$ distributions (Cossu et al. Reference Cossu, Pujals and Depardon2009; Fan et al. Reference Fan, Kozul, Li and Sandberg2024). Pickering et al. (Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021) demonstrated that the classical eddy-viscosity models (e.g. RANS) behave identically to an optimal eddy-viscosity model in aligning the most energetic mode. This agreement also underscores that the predicted most energetic response and scales derived from RANS-based $ \mu _{t}$ , despite inherent modelling uncertainties, retain physical fidelity and provide reliable predictions. Consequently, the subsequent parametric investigations in this study employ RANS-based flows, balancing computational efficiency with demonstrated predictive fidelity.

5.2. Geometric parameters

Two geometric parameters are considered: the total turning angle $ \varphi$ , and the curvature radius $ r$ . It should be noted that the frequency is set to $ St=0.036$ , which is a representative frequency for Görtler vortices.

Figure 17.

( $ a$ ) Normalised optimal gains at different turning angles as functions of spanwise wavenumbers $ \beta \delta$ . ( $ b$ ) Distributions of integrated Chu energy density along the wall-normal direction at different turning angles at their respectively optimal wavenumber. The inset in ( $ a$ ) represents the real part of $ w^{\prime}$ for the optimal response in the case $ \varphi =30^{\circ }$ , while the arrows indicate local maxima for different cases. The black dash-dotted line in ( $ b$ ) indicates the beginning of the concave surfaces, the black dashed line represents the predicted increasing trend by LSA in figure 7( $ a$ ), and the open circles denote the ends of the concave surfaces.

Figure 17( $ a$ ) shows the distributions of the optimal gains at different turning angles for a constant curvature radius $ r=50\delta $ . Three peaks are marked in figure 17( $ a$ ), corresponding to $ \varphi =0^{\circ }$ , $ \varphi =5^{\circ }$ and $ \varphi \geqslant 10^{\circ }$ . No centrifugal effects are excited at $ \varphi =0^{\circ }$ because it belongs to the flat plate. Consequently, the peak wavenumber $ \beta \delta =1.86$ corresponds to coherent structures in zero pressure gradient cases, as reported by Cossu et al. (Reference Cossu, Pujals and Depardon2009), Alizard et al. (Reference Alizard, Pirozzoli, Bernardini and Grasso2015) and Abreu et al. (Reference Abreu, Cavalieri, Schlatter, Vinuesa and Henningson2020). We do not discuss coherent structures in flat plate cases further, as there are no centrifugal effects. When the total turning angle increases to $ \varphi =5^{\circ }$ , weak centrifugal effects cause the optimal wavenumber to shift to $ \beta \delta =2.59$ . As the total turning angle increases to $ \varphi =10^{\circ }$ , the optimal wavenumber shifts again due to the elongated region affected by the centrifugal effects. Further increasing the turning angle no longer changes the optimal wavenumber, which indicates that the optimal wavenumber is unrelated to the turning angle when $ \varphi \geqslant 10^{\circ }$ . This observation is consistent with our modal analysis results, where the spanwise scale changes slightly at the middle and end of the concave surface, as shown in figures 11 and 15. There may be a critical angle for a specific concave radius $ r$ beyond which the optimal wavenumber is independent of the turning angle. The critical value is approximately $ \varphi =10^{\circ }$ for $r = 50 \delta$ . Additionally, changes in the total turning angle at the same concave radius do not alter the coherent structures qualitatively, as shown in the inset of figure 17( $ a$ ) in the case $ \varphi =30^{\circ }$ . Notably, no flow separation occurs in these cases.

Figure 17( $ b$ ) shows the distributions of the integrated Chu energy density along the wall-normal direction for different total turning angles at their respective optimal wavenumbers. A direct comparison of the Chu energy across different cases is reasonable because the Chu energy is normalised by the energy of the forcing. The Chu energy densities exhibit a similar exponential increase in the concave region, which agrees well with the spatial growth rate predicted by LSA. This common increasing trend is due to the constant streamline concave curvature $ K$ in the concave region, while different turning angles have different Chu energy levels at the ends of the concave surfaces.

Figure 18.

( $ a$ ) Optimal wavenumbers and ( $ b$ ) corresponding distributions of integrated Chu energy density at different concave curvatures. The insets in ( $ a$ ) represent the real parts of $ w^{\prime}$ of the optimal responses, in the cases $ K=0.014$ and $ K=0.067$ , and the red dashed line indicates $ \beta \delta =3.2$ . The black dashed line in ( $ b$ ) marks the beginning of the concave surface, while the open circles in ( $ b$ ) represent the ends of the concave surfaces.

The second geometric parameter is the curvature radius $ r$ , which corresponds to the concave curvature $ K$ . Figure 18( $ a$ ) illustrates the optimal wavenumbers at different concave curvatures for a constant total turning angle $ \varphi =20^{\circ }$ . The optimal wavenumber rises sharply at small $ K$ and more gradually at larger values, with only a minor increase between $ K=0.083$ and $ K=0.1$ . It is predicted that further increasing $ K$ will cause the optimal wavenumber to converge towards a specific constant value, approximately $ \beta \delta =3.2$ , as indicated by the red dashed line in the figure. The coherent structures in the cases $ K=0.014$ and $ K=0.067$ shown in the insets of figure 18( $ a$ ) resemble those at $ K =0.02$ in figure 5, which implies that changes in concave curvature do not affect the coherent structures qualitatively.

Figure 18( $ b$ ) presents the integrated Chu energy density along the wall-normal direction at different concave curvatures at their respective optimal wavenumbers. The exponentially increasing trends in the concave region are evident for cases at different curvatures. Unlike the consistent trend observed in figure 17( $ b$ ), the trends differ for various curvatures. Specifically, a large curvature results in a rapid increase in Chu energy density. The rapid increase of fluctuation energy in a large concave curvature region aligns with previous experimental findings (Jayaram et al. Reference Jayaram, Taylor and Smits1987). Although the region influenced by the centrifugal effects is short at a large curvature, the Chu energy density still reaches a high level.

An infinite curvature $ K\rightarrow \infty$ represents a $ 20^{\circ }$ compression ramp configuration. However, achieving $ K\rightarrow \infty$ for a streamline in real flows is impossible because an adverse pressure gradient will cause a separation bubble that alters the curvature. Figure 19( $ a$ ) shows the concave curvature distributions of three streamlines passing through $(x_0,0.01\delta)$ , $(x_0,0.05\delta)$ and $(x_0,0.15\delta)$ in a $ 20^{\circ }$ compression ramp case. The concave curvatures are at large values near the separation and reattachment points, but remain small close to the corner. The two pronounced curvature peaks near these points are significantly larger than the critical value 0.01 (Hoffmann et al. Reference Hoffmann, Muck and Bradshaw1985), making the optimal wavenumber approach the upper boundary $ \beta \delta =3.2$ under the individual influence of the centrifugal effects. However, in a dynamic system of turbulent flows with a separation bubble, bubble dynamics and the centrifugal effects coexist. The bubble dynamics may dominate the centrifugal effects on the main amplification mechanism in the dynamic system predicted by resolvent analysis. As mentioned by Hao (Reference Hao2023), in a $ 25^{\circ }$ compression ramp case, the local maxima of the optimal gain are related to the bubble dynamics, while the Görtler instability is insignificant.

Figure 19.

( $ a$ ) Distributions of concave curvatures of three streamlines passing through $(x_0,0.01\delta)$ , $(x_0,0.05\delta)$ and $(x_0,0.15\delta)$ in a $ 20^{\circ }$ compression ramp case. ( $ b$ ) Normalised optimal gains as functions of spanwise wavenumbers $ \beta \delta$ at $ St=0.036$ for different geometries. The black lines in ( $ a$ ) denote the mean separation and reattachment points. Local maxima at $ \beta \delta =3.3$ are indicated by the black dashed line in ( $ b$ ).

To validate the predicted upper boundary value and illustrate the large concave curvature effects, the three streamlines in figure 19( $ a$ ) are established as walls of new geometries to isolate the separation bubble. For comparison, an additional geometry with a straight line connecting the separation and reattachment points is also considered. Figure 19( $ b$ ) compares the optimal gains obtained using the four new geometries. As expected, all four cases share an analogous trend and a common peak $ \beta \delta =3.3$ , i.e. $ \lambda _{z}=1.9\delta$ , which is close to the $ 2\delta$ Görtler structures reported by experiments (Schuelein & Trofimov Reference Schuelein and Trofimov2011; Zhuang et al. Reference Zhuang, Tan, Liu, Zhang and Ling2017) and numerical simulations (Loginov, Adams & Zheltovodov Reference Loginov, Adams and Zheltovodov2006; Priebe et al. Reference Priebe, Tu, Rowley and Martín2016) in compression ramp cases. Furthermore, the Görtler structures are also found in impinging shock interactions (Pasquariello, Hickel & Adams Reference Pasquariello, Hickel and Adams2017; Li et al. Reference Li, Zhang, Yu, Lin, Tan and Sun2022) where separation bubbles also exist. Obviously, the coupled effects of separation bubbles with the Görtler instability influence the Görtler structures slightly. Therefore, the upper boundary value provides some valuable insights into characterising the spanwise scales of Görtler vortices with separation bubbles.

One might argue that the shock structures of the four cases are significantly different from those with the separation bubble in the $ 20^{\circ }$ compression ramp case. However, an equal strength of the concave curvature is the essential parameter for separately investigating the Görtler vortices in compression ramp cases.

5.3. Freestream parameters

Additional cases are considered to investigate the influence of freestream parameters on the spanwise scales of Görtler vortices, including the unit Reynolds number $ {Re}_{\infty }$ , Mach number $ {M}_{\infty }$ , and temperature ratio $ {T}_{w}/{T}_{\infty }$ . We isolate the influence of one parameter by keeping the other two parameters unchanged. The cases with the original freestream parameters are referred to as the baseline cases.

Figure 20.

Normalised optimal gains as functions of spanwise wavenumbers $ \beta \delta$ at different ( $ a$ ) Mach number, ( $ b$ ) temperature ratio, and ( $ c$ ) Reynolds number, at $ r=50\delta$ and $ \varphi =20^{\circ }$ . The red arrows in (a,b) indicate increasing trends, and the black dashed line in ( $ c$ ) means a local maximum.

Figure 21.

Preferential wavenumbers as a function of concave curvature for a constant turning angle $ \varphi =20^{\circ }$ .

Optimal gains as a function of the scaled wavenumber $ \beta \delta$ at different freestream parameters are shown in figure 20. The optimal wavenumber increases with an increase in $ {M}_{\infty }$ . A similar increasing trend is shown in figure 20( $ b$ ) when decreasing the temperature ratio $ {T}_{w}/{T}_{\infty }$ , although the change is slight. Figure 20( $ c$ ) depicts the optimal gains at different freestream Reynolds numbers $ Re_{\infty }$ . The optimal wavenumber $ \beta \delta =2.73$ is shared by different Reynolds number cases, implying that the optimal wavenumber predicted by resolvent analysis is independent of the freestream Reynolds number, when normalised using the boundary-layer thickness.

Figure 21 shows the optimal wavenumbers as a function of concave curvature under different flow conditions. Trends similar to the baseline cases are observed at different flow parameters, with a rapid rise when $ K$ is small, and a slow rise when $ K$ is large. The upper boundary corresponding to $ K\rightarrow \infty$ can also be determined for various freestream parameters. As discussed in § 5.2, the upper boundary value can be regarded as an approximate wavenumber of the Görtler vortices in cases with a separation bubble. Hence the Mach number effect can be verified qualitatively by different supersonic or hypersonic compression ramp cases. At $ {M}_{\infty } =2.95$ , the structures of $ 2\delta$ were observed from time-averaged skin friction distribution by Loginov et al. (Reference Loginov, Adams and Zheltovodov2006). When $ {M}_{\infty } =6$ , the Görtler structures downstream of the separation bubble are changed to $ 1.2\delta$ , as shown by Guo et al. (Reference Guo, Zhang, Zhu and Li2024). The decrease in the wavelength of Görtler vortices is consistent with our findings. The decreasing trend may be related to compressibility effects, which will be investigated in a future study. Furthermore, the Reynolds number independent relation is verified, as the optimal wavenumbers overlap for the two Reynolds number cases at different curvatures, in agreement with experimental findings by Schuelein & Trofimov (Reference Schuelein and Trofimov2011) in turbulent separated flows.

Our results offer quantitative insights into predicting the scale of Görtler vortices. Future efforts should focus on exploring the coupling effects between separation bubbles and Görtler instability in compression ramps, as well as impinging shock interactions in turbulent flows.