Characterising Görtler vortices in supersonic turbulent flows over concave surfaces

Abstract

Görtler vortices induced by concave curvature in supersonic turbulent flows are investigated using resolvent analysis and large-eddy simulations at Mach 2.95 and Reynolds number $ Re_{\delta }=63\,500$ based on the boundary-layer thickness $ \delta$. Resolvent analysis reveals that the most amplified coherent structures manifest as streamwise counter-rotating vortices with optimal spanwise wavelength $ 2.4\delta$ at cut-off frequency $f\delta /{u}_{\infty } =0.036$, where $ {u}_{\infty }$ is the freestream velocity. The leading spectral proper orthogonal decomposition modes with spanwise wavelength approximately $ 2\delta$ align well with the predicted coherent structures from resolvent analysis at $f\delta /{u}_{\infty } =0.036$. These predicted and extracted coherent structures are identified as Görtler vortices, driven by the Görtler instability. The preferential spanwise scale of the Görtler vortices is further examined under varying geometric and freestream parameters. The optimal spanwise wavelength is insensitive to the total turning angle beyond a critical value, but sensitive to the concave curvature $ K$ at the same turning angle. A limit spanwise wavelength $ 1.96\delta$, corresponding to an infinite concave curvature as $ K\rightarrow \infty$, is identified and validated. Increasing the freestream Mach number or decreasing the ratio of wall temperature to freestream temperature reduces the optimal wavelength normalised by $ \delta$, while variations in freestream Reynolds number have negligible impact. Additionally, a modified definition of the turbulent Görtler number $ G_{T}$ based on the peak eddy viscosity in boundary layers is proposed and employed to assess the occurrence of Görtler instability.

1. Introduction

Compressible turbulent boundary layers (TBLs) over concave surfaces are crucial in configurations such as aerofoils and scramjet engines. The combined effects of concave streamline curvature, adverse pressure gradients and bulk compression in flows over concave surfaces destabilise TBLs to promote turbulent mixing, as highlighted by Bradshaw & Young (1973) and Bradshaw (1974). According to Floryan (1991), the excited centrifugal effects affect TBLs from three aspects: (i) directly influencing turbulent structures by the wall curvature effect, (ii) generating secondary flows known as Görtler vortices, and (iii) impacting turbulent structures via Görtler vortices. The second effect is also referred to as Görtler instability. However, the behaviours and characteristics of the induced Görtler vortices in turbulent flows remain inadequately understood compared to those in laminar flows.

Early studies focused on experiments to explore the impact of concave curvature $ K=\delta /r$ on turbulence, where $ \delta$ is the boundary-layer thickness, and $ r$ is the curvature radius. Hoffmann, Muck & Bradshaw (1985) experimentally investigated the impact of curvature on incompressible TBLs. They concluded that the wall curvature directly dominates turbulent structures when the curvature is not strong enough, while the Görtler vortices are more influential at a moderate curvature $ K=\delta /r\geqslant 0.01$ . The critical value suggests that strong curvature is necessary to form Görtler vortices due to fuller TBLs (Floryan 1991). In supersonic flows, Jayaram, Taylor & Smits (1987) compared TBLs over concave surfaces with curvatures 0.02 and 0.1 to those over a compression ramp ( $ K\rightarrow \infty$ ) with the same turning angle $ 8^\circ$ . The study revealed that turbulent fluctuations respond rapidly in regions with large curvature to alter turbulent structures. The authors also reported a dip below the logarithmic law, but did not detect longitudinal roll cells in their experiments. Later, Barlow & Johnston (1988a , b ) observed that longitudinal roll cells wander side to side without preferred spanwise locations when the freestream is relatively free of spanwise non-uniformities, indicating that these roll cells are unsteady. However, a stationary pattern of large-scale roll cells was found when introducing steady weak disturbances using small vortex generators. These stable structures were observed in time-averaged visualisations with the spanwise wavelength approximately 1–2 times the boundary-layer thickness. A similar steady pattern with the help of upstream vortex generators was also discussed by Schuelein & Trofimov (2011) in supersonic compression ramp cases. Different patterns of the large-scale motions highlight the influence of upstream disturbances, as summarised by Floryan (1991) and Priebe et al. (2016). Donovan, Spina & Smits (1994) investigated the response of the large-scale motions to a short concave region but with a strong concave curvature $ K = 0.08$ . They reported that the centrifugal effects amplify turbulence levels and skin friction compared to flat plates with similar pressure distributions. Later, Wang & Wang (2016) analysed the break-up process in the concave region from large-scale structures to smaller ones with high-resolution flow visualisation, finding that the vortices tend to break up faster under more significant concave curvature.

In recent years, high-fidelity numerical simulations have been employed to investigate turbulent structures over concave surfaces. Tong et al. (2017) performed direct numerical simulations (DNS) to study supersonic TBLs over a curved ramp with $ K=0.02$ and turning angle $ \varphi =24^{\circ }$ at Mach 2.95. They discovered a mode at $ f\delta /{u}_{\infty } =0.017$ with regular spatial structures distributed in the spanwise direction using dynamic mode decomposition, where $ {u}_{\infty }$ is the freestream velocity. These structures resemble the spanwise structures of Görtler vortices in laminar flows. Sun, Sandham & Hu (2019) compared TBLs over concave surfaces with two different radii to those over a flat plate, and reported prominent low-frequency fluctuations (within the range $ f\delta /{u}_{\infty } =0.2{-}0.4$ ) on concave surfaces due to large-scale structures. They noted that Görtler instability primarily affects the outer layer, forming Görtler-like vortices, while the inner layer remains largely unaffected. These vortices facilitate rapid momentum exchange between the inner and outer parts of TBLs, enhancing turbulence intensity and Reynolds stresses. Wang et al. (2019) also observed similar momentum exchange processes driven by large-scale roll cells, achieved through ejection and sweep events that transport high-momentum flow towards the wall, and low-momentum flow towards the outer edge. They suggested that these large-scale roll cells originate from very-large-scale motions in the absence of upstream artificial disturbances. Wu, Liang & Zhao (2019) numerically investigated the small curvature case studied by Jayaram et al. (1987) at $ K=0.02$ and $ \varphi =8^{\circ }$ , but with a lower reference Reynolds number. They found that the outer peak of the energy spectra falls in the lower-wake region due to the turbulence amplification effect. However, they demonstrated that log-region superstructures primarily contribute to outer–inner interactions in the concave region, i.e. large/small-scale interactions, enhancing turbulence modulation strength significantly. In a recent study, Tong et al. (2023) carried out DNS to characterise the effect of concave curvature on the wall shear stress and heat flux. They emphasised the leading role of the outer large-scale structures in contributing to the mean wall shear stress and wall heat flux.

Previous studies have primarily examined the effects of concave curvature on turbulent statistical properties, emphasising how wall curvature and Görtler vortices significantly alter turbulent structures. However, the Görtler vortices themselves in turbulent flows, particularly their spatial–temporal structures and scales, remain poorly understood. These vortices, characterised by counter-rotating roll cells in laminar flows, are challenging to visualise in turbulent flows due to the presence of complex small-scale structures (Wang et al. 2019). Additionally, while many studies have shown that the spanwise wavelength of Görtler vortices is typically approximately $ 2\delta$ , Guo et al. (2024) also reported structures with spanwise wavelength approximately $ 1.2\delta$ . Therefore, the spanwise scales of Görtler vortices also need to be thoroughly examined.

In recent years, resolvent analysis (McKeon & Sharma 2010) has become an attractive tool for studying coherent structures in turbulent flows. This method identifies the most amplified response for a given input forcing, and has been successfully applied to predict coherent structures in turbulent jets (Schmidt et al. 2018; Lesshafft et al. 2019; Maia et al. 2024) and channel flows (McKeon 2017; Bae, Dawson & McKeon 2020; Fan et al. 2024), as well as to investigate the low-frequency dynamics of turbulent separation bubbles (Hao 2023; Cura et al. 2024).

The present study aims to investigate the features of Görtler vortices in supersonic turbulent flows using resolvent analysis and large-eddy simulations (LES). We will focus on the spatial–temporal structures and spanwise scales of these vortices. Resolvent analysis will be used to predict the most amplified coherent structures on concave surfaces. These predictions will be compared with structures extracted by proper orthogonal decomposition (POD) and spectral proper orthogonal decomposition (SPOD) based on the LES dataset. We expect the predicted structures to be Görtler vortices, with Görtler instability as the dominant amplification mechanism. Finally, we will explore how geometric and freestream parameters influence the preferential spanwise scale of Görtler vortices and predict their occurrence.

2. Geometry and freestream conditions

As sketched in figure 1, the flow configuration in this study involves a concave surface with total turning angle $ \varphi =20^{\circ }$ and curvature radius $ r=50 \delta$ ( $ K=0.02$ ). A Cartesian coordinate system is established with the origin at the beginning of the concave surface. The computational domains for the Reynolds-averaged Navier–Stokes (RANS) simulations and LES are also marked in the figure. The freestream parameters are taken from Zheltovodov et al. (1990), with Mach number $ {M}_{\infty }=2.95$ , density $ {\rho }_{\infty }=0.314$ kg m $^{-3}$ , statistic temperature $ {T}_{\infty }=108$ K, and Reynolds number $ {Re}_{\delta }=63\,500$ with $ \delta = 2.27$ mm. The boundary-layer thickness $ \delta$ is defined as $ 99\, \%$ of the freestream velocity at the reference station $ x_{0}$ (see figure 1).

Figure 1. Computational domains and boundary conditions: $ r$ is the curvature radius, $ \varphi$ is the turning angle, $ x_{0}$ is at $ (-1.76\delta ,0)$ , the reference station, and $ x_{in}$ is at $(-22.5\delta ,0)$ , the station of the chosen one-dimensional velocity profile from the RANS base flow. The RANS domain and LES domain refer to the computation domains for RANS simulations and LES, respectively. Boundary conditions for RANS and LES are also marked.

3. Theoretical and numerical approach

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 2009).

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 1992). An in-house finite-volume solver called PHAROS (Hao, Wang & Lee 2016; Hao & Wen 2020) is used to obtain RANS base flows. The modified Steger–Warming scheme (MacCormack 2014) 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 (1998) 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 2010) combining the AUSM $ ^{+}$ -up scheme (Liou 2006) with the fifth-order weighted essentially non-oscillatory reconstruction (Jiang & Shu 1996) is used to solve the inviscid fluxes, switched by the Jameson sensor (Jameson, Schmidt & Turkel 1981). 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 2009; Ceci et al. 2022) 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 2012). 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 1990) and investigated using resolvent analysis (Sipp & Marquet 2013). 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. 2019). 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 1965) is applied to the response and the forcing to evaluate their respective energies. As discussed by Bugeat et al. (2019), (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 2018). 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 2012). Finally, the implicitly restarted Arnoldi method implemented in ARPACK (Lehoucq, Sorensen & Yang 1996), interfaced with the MUMPS (Amestoy et al. 2001) 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. 2023), 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. 2024). 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. (2019). 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. (2024).

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 2014) 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 & 2010) 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 1967; Towne, Schmidt & Colonius 2018; Schmidt et al. 2018) 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. 2018). 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 1965) 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. 2018).

4. Spatial–temporal structures

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 2009). 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 2001). 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 2001), 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 2001). 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. 2020; Hao et al. 2023). 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. 1985). 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 1996; Tumin & Reshotko 2003). 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. 2019). 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. 2019). Similar distributions in $ k$ and Reynolds stress components were also discussed by Wu et al. (2019) and Sun et al. (2019).

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 1965). 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. 2018; Mendez, Balabane & Buchlin 2019). 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 1967) 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 1965). A value $ \gamma = 0$ indicates orthogonality between the two modes, while $ \gamma = 1$ signifies perfect alignment (Abreu et al. 2020; Fan et al. 2024). 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. (2024) 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. (2018), Schmidt et al. (2018) and Abreu et al. (2020). When forcings are unit-variance white noise, SPOD is equivalent to resolvent analysis. In real turbulent flows, biased and nonlinear forcings (Morra et al. 2021) 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. 2020). 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 2024). 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. 2021). 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 2021; Yeung & Schmidt 2024). 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. (2017). 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 1988a , b ; Schuelein & Trofimov 2011).

5. Scaling of Görtler vortices

The alignment between the leading SPOD modes obtained from the LES dataset and resolvent modes based on the LES base flow confirms the validity of resolvent analysis for capturing Görtler vortices dynamics. While resolvent analysis with LES base flows achieves high fidelity, a systematic parametric investigation of the spanwise scale sensitivity to geometric and freestream parameters – across multiple configurations – would be computationally prohibitive if relying solely on LES-based resolvent analysis. To balance accuracy and practicality, we adopt RANS-based resolvent analysis for the parametric study in this section. To validate this approach, § 5.1 first compares RANS-based resolvent analysis results against the LES-based methodology, demonstrating quantitative agreement in mode structures and optimal gains. Subsequent subsections then leverage RANS-based analysis to examine parameter dependencies, with the validated methodology ensuring robustness of the conclusions.

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 (2009). 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. 2009; Fan et al. 2024). Pickering et al. (2021) 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. (2009), Alizard et al. (2015) and Abreu et al. (2020). 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. 1987). 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. 1985), 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 (2023), 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 2011; Zhuang et al. 2017) and numerical simulations (Loginov, Adams & Zheltovodov 2006; Priebe et al. 2016) in compression ramp cases. Furthermore, the Görtler structures are also found in impinging shock interactions (Pasquariello, Hickel & Adams 2017; Li et al. 2022) 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. (2006). 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. (2024). 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 (2011) 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.

6. Occurrence of Görtler instability

The results from resolvent analysis and LSA are examined in § 4.1. The LSA provides accurate predictions for both the spatial growth rate and the spanwise wavelength of the Görtler mode, consistent with values predicted by resolvent analysis. Consequently, LSA can be employed to estimate the spatial growth rate of the Görtler mode, thereby assessing the occurrence of Görtler vortices.

Two parameters are usually used to evaluate Görtler instability in turbulent flows: the concave curvature $ K$ , and the Görtler number $ {G}_{T}$ . The definition of $ {G}_{T}$ in laminar flows is given by

(6.1) \begin{equation} {G}_{T}=\frac {{\rho }_{\infty } {u}_{\infty }}{{\mu }_{\infty } }\theta \sqrt {\frac {\theta }{r}}, \end{equation}

where $ \theta$ represents the boundary-layer momentum thickness. Tani (1962) suggested using the eddy viscosity $ \mu _{t}$ to replace the molecular viscosity $ {\mu }_{\infty }$ for defining $ {G}_{T}$ in turbulent flows. He followed a simple assumption of a constant eddy viscosity $ \mu _{t}=0.018 {\rho }_{\infty } {u}_{\infty }\delta ^{*}$ for outer layers of incompressible TBLs (Clauser 1956), where $ {\delta }^{*}$ is the boundary-layer displacement thickness. This leads to the following expression for the turbulent Görtler number:

(6.2) \begin{equation} {G}_{T}=\frac {\theta }{0.018\delta ^{*} } \sqrt {\frac {\theta }{r}} . \end{equation}

Smits & Dussauge (2006) extended the turbulent $ {G}_{T}$ to compressible flows while using the same assumption of a constant eddy viscosity. However, the assumption of a constant $ \mu _{t}$ in outer layers of TBLs may not be appropriate. One can refer to Abe & Antonia (2019) and figure 16( $ a$ ) for distributions of $ \mu _{t}$ in incompressible and compressible flows, respectively. Furthermore, the empirical formula overestimates $ \mu _{t}$ in the outer layer of TBLs. For instance, the formula gives $ {\mu }_{t}/{\mu }_{\infty }=365$ , which is significantly higher than the peak values $ {\mu }_{t}/{\mu }_{\infty }=125$ and 86 at the reference station $ x_{0}$ from RANS and LES. Such an overestimation causes an underestimation of the turbulent $ {G}_{T}$ . Different modelling approaches for $ \mu _{t}$ based on DNS/LES databases and RANS models of zero pressure gradient flat plates lead to slight variations in the peak $ \mu _{t, peak}$ values in TBLs. To investigate the critical $ G_{T}$ under different flow conditions, we choose the peak value $ \mu _{t, peak}$ at the reference station $ x_{0}$ from the RANS results (see figure 16 $ a$ ) as the reference value to replace the $ {\mu }_{\infty }$ in calculating the turbulent $ {G}_{T}$ , given as

(6.3) \begin{equation} {G}_{T}=\frac {{\rho }_{\infty } {u}_{\infty }}{{\mu }_{t, peak } }\theta \sqrt {\frac {\theta }{r}}. \end{equation}

As we mentioned before, the critical value $ K=0.01$ for inducing streamwise vortices was proposed by Hoffmann et al. (1985), and used by Loginov et al. (2006) and Priebe et al. (2016) in turbulent separated flows. Smits & Dussauge (2006) also suggested that the critical value $ K$ is approximately 0.03 in Mach 3 flows, and changes as Mach number changes. However, our LES results clearly show that Görtler vortices begin to manifest at $ K=0.02$ , indicating that the value 0.03 may be somewhat overestimated.

Moreover, the critical $ {G}_{T}$ for turbulent flows has not been established (Floryan 1991). Many researchers (Loginov et al. 2006; Priebe et al. 2016; Tong et al. 2017) have used the critical value 0.6 in laminar flows (Görtler 1954) as a criterion for turbulent flows (calculated from (6.2)). However, the value would increase to a range 1.8–2.4 if calculated using the correct $ {\mu }_{t}$ . Given that the growth rate predicted by LSA agrees well with the resolvent analysis result, LSA is used in this study to determine the critical $ {G}_{T}$ based on (6.3).

Figure 22. ( $ a$ ) The maximum spatial growth rate $ -\alpha _{i} \delta$ at station $ \varphi = 10^{\circ }$ predicted by LSA at $ St=0.036$ . ( $ b$ ) The optimal wavenumber $ \beta \theta$ from resolvent analysis, along with the most unstable wavenumber $ \beta \theta$ from LSA of the baseline cases as a function of Görtler number $ G_{T}$ . The Görtler numbers are extracted from streamlines passing through $(x_0,0.5\delta)$ . The black dashed line in ( $ a$ ) indicates zero growth rate. The black line in ( $ b$ ) indicates $ \Lambda =215$ , where $ \Lambda =( {{\rho }_{\infty } {u}_{\infty }}/{{\mu }_{t, peak } }){\lambda }_{z}\sqrt {{{\lambda }_{z}}/{r}}$ is a dimensionless wavelength parameter (Saric et al. 1994).

Figure 22( $ a$ ) presents the maximum spatial growth rates at station $ \varphi = 10^{\circ }$ predicted by LSA, with the peak $ {G}_{T}$ of streamlines passing through $(x_0,0.5\delta)$ for the baseline cases across varying concave curvatures. It can be observed that as $ {G}_{T}$ progressively diminishes, the maximum growth rate $ -\alpha _{i} \delta$ exhibits a concomitant decline. A maximum spatial growth rate $ -\alpha _{i} \delta \lt 0$ signifies that the Görtler mode remains unexcited, resulting in the absence of Görtler vortices. A zero growth rate is attained at $ {G}_{T}= 0.3$ , corresponding to $ K= 0.00143$ . The predicted $ {G}_{T}=0.3$ is close to the laminar critical value 0.6 (Görtler 1954), but the curvature $ K=0.00143$ is significantly lower than the experimentally obtained empirical value 0.01 for turbulent flows (Hoffmann et al. 1985). The critical $ G_{T}$ becomes 0.45 when using $ \mu _{t, peak}$ from LES. Both values are much lower than the previous empirical range 1.8–2.4 (calculated using (6.3)). It should be noted that different modelling methods for $ \mu _{t}$ may lead to a slight change in the critical $ G_{T}$ values, although the trends remain similar. Essentially, the critical value 0.6 in laminar flows (or 0.4638 in Floryan & Saric 1982) is an empirical guideline at the incompressible limit. However, El-Hady & Verma (1984) found that the critical $ G_{T}$ increases with rising $ M_{\infty }$ , which indicates that this guideline may not be suitable for all laminar flow scenarios, such as varying Mach numbers and wall temperatures. There may not be a solid criterion for compressible laminar and turbulent flows. Thus we do not aim to specify a precise critical $ G_{T}$ , but rather provide an approximate range 0.3–0.45 to offer meaningful guidance on the occurrence of Görtler vortices in compressible turbulent flows. Future work will incorporate multi-fidelity $ {\mu }_{t}$ from DNS/LES datasets to refine this range, but our current results provide a conservative estimate for design purposes.

It is worth noting that the critical values obtained, whether for laminar or turbulent flows, are only of reference significance. Only when $ G_{T}$ is much larger than the critical value will the Görtler vortices manifest obviously in flows. In practice, determining the critical $ G_{T}$ experimentally poses significant challenges due to the low energy growth rates, which necessitates a long concave surface to allow Görtler vortices to develop observable amplitudes. However, the finite length of models artificially truncates this growth process, leading to overestimated critical $ G_{T}$ values – often exceeding theoretical predictions by orders of magnitude. Consequently, Görtler vortices may remain undetected even in configurations where $ G_{T}$ far surpasses the theoretical threshold. The determined $ G_{T}$ and $ K$ by LSA represent theoretical lower bounds.

Figure 22( $b$ ) compares the optimal wavenumber from resolvent analysis with the most unstable wavenumber obtained from LSA as a function of $ {G}_{T}$ for the baseline cases. When $ {G}_{T}$ is small, i.e. the concave curvature is small, the wavenumbers from both methods align closely with each other. Notably, both analyses yield wavenumbers that fall along the line $ \Lambda =215$ when $ \beta \theta \lt 0.18$ , which is consistent with the constant- $ \Lambda$ assumption in laminar flows (Tani 1962; Saric et al. 1994). However, a discrepancy occurs when $ {G}_{T}\gt 1.2$ . In this case, the most unstable wavenumber from LSA remains on the line $ \Lambda =215$ , while the optimal wavenumber from resolvent analysis gradually converges to a constant value as $ {G}_{T}$ increases, in agreement with the upper boundary value in figure 18. This discrepancy is attributed to non-parallelism and non-modal effects (Hanifi et al. 1996; Tumin & Reshotko 2003). Additionally, figure 22( $b$ ) also provides an approximation for the spanwise wavenumbers of Görtler vortices within the Görtler number range 0.3–1.2.

7. Conclusions

The spatial–temporal features of Görtler vortices in supersonic turbulent flows over concave surfaces are investigated using resolvent analysis and LES. Resolvent analysis predicts an optimal response with the spanwise wavelength $ \lambda _{z} \approx 2.4\delta$ and the cut-off frequency $St=0.036$ , based on an LES base flow at Mach 2.95 and $ Re_{\delta }=63\,500$ over a concave surface with a total turning angle $ \varphi =20^{\circ }$ and curvature radius $ r=50 \delta$ . The optimal resolvent mode manifests as streamwise counter-rotating vortices in the concave region, identified as the Görtler vortices, which are consistent with the spatial structures in laminar flows (Görtler 1954; Floryan 1991). The dominant amplification mechanism is the Görtler instability. The integrated Chu energy density increases exponentially in the concave region, with the growth rate closely matching LSA results. Furthermore, the optimal wavelength and mode shapes predicted by LSA also align well with those obtained by resolvent analysis.

These predictions are validated against the same LES dataset. The centrifugal effects promote momentum exchange and enhance turbulence intensity, with a second peak of turbulent kinetic energy occurring at $y_{n}/\delta = 0.26$ . The POD analysis is performed on $y$ – $z$ slices at stations $\varphi = 10^{\circ }$ and $20^{\circ }$ . The most energetic spatial structures at the two stations are similar streamwise counter-rotating vortices with spanwise wavelength approximately $ 2\delta$ . Additionally, coherent structures in the spanwise mid-plane and several $y$ – $z$ planes are extracted using SPOD. The leading SPOD mode of the spanwise mid-plane peaks at $ St = 0.036$ , resembling the optimal output from resolvent analysis at the same frequency. The highest projection coefficient between the leading SPOD mode and the optimal resolvent response is approximately 0.94. The leading SPOD modes at three streamwise stations $\varphi = 10^{\circ }$ , $\varphi = 15^{\circ }$ and $\varphi = 20^{\circ }$ also show similar spanwise structures around $2\delta$ , consistent with the reconstructed structures using the optimal resolvent mode. The good agreement between the optimal SPOD modes and the resolvent mode implies that resolvent analysis accurately captures the Görtler vortices.

Further investigations are conducted to reveal how the preferential spanwise wavelength of Görtler vortices scales with varying geometric and freestream parameters. The optimal wavelengths are insensitive to the total turning angle beyond a critical value. However, they change as the concave curvature varies at the same turning angle. A limit wavelength $ \lambda _{z}=1.96\delta$ ( $\beta \delta =3.2$ ), corresponding to infinite curvature $ K\rightarrow \infty$ , is identified and validated using a $ 20^{\circ }$ compression ramp case by isolating the separation bubble. This value provides valuable insights into the Görtler instability in turbulent separated flows. The limit wavelength decreases with an increasing Mach number $ {M}_{\infty }$ and a decreasing temperature ratio $ {T}_{w}/{T}_{\infty }$ , remaining unchanged with varying Reynolds number $ Re_{\infty }$ when normalised using the boundary-layer thickness $ \delta$ . Additionally, a modified definition of the turbulent Görtler number $ G_{T}$ is proposed using the peak eddy viscosity in TBLs. The critical Görtler number $ G_{T}=0.3{-}0.45$ and curvature $ K=0.00143$ are determined by LSA and serve as conservative estimates for evaluating the occurrence of Görtler vortices in the baseline cases.

The optimal forcings responses may provide guidance for flow control strategies. For instance, these optimal forcings can be used through plasma actuators (Ribeiro & Taira 2024) to suppress or enhance the strengths of Görtler vortices. Future studies will discuss control strategies in detail.