On the low-frequency unsteadiness in shock wave–turbulent boundary layer interactions

Abstract The shock wave–turbulent boundary layer interaction over a compression corner is studied using global stability analysis (GSA) and resolvent analysis based on a separation of scales between the low-frequency, large-scale motions and the turbulent fluctuations. The GSA identifies a leading stationary mode, which becomes globally unstable as the ramp angle is beyond a critical value. For globally stable flows, the resolvent analysis captures two-dimensional and three-dimensional local maxima in optimal gain, both of which are due to modal resonance between the forcing and the leading global mode. Notably, the frequency-premultiplied optimal gain associated with two-dimensional disturbances peaks at a low frequency. For different interaction strengths, the peak frequencies collapse onto a universal value of 0.015 when non-dimensionalized using the length of the separation region and the free-stream velocity. A numerical simulation perturbed with the corresponding optimal forcing reveals that the response is in the form of a back-and-forth shock motion.


Introduction
A shock wave-turbulent boundary layer interaction (STBLI) is a fundamental problem of fluid mechanics, which is also of practical importance in the aerodynamics and propulsion of high-speed flight vehicles (Clemens & Narayanaswamy 2014).A remarkable feature of supersonic STBLIs is a low-frequency, large-scale motion of the flow structure, the source of which is still an open issue.
In addition to the well-known low-frequency back-and-forth motion of the separation shock, the separation bubble exhibits large-scale three-dimensionality. Settles, Fitzpatrick Magidov 2007; Crouch et al. 2009;He & Timme 2021).In contrast, a noise amplifier feeds on continuous input of external disturbances.Resolvent analysis is usually used to study the dynamics of noise amplifiers, which has also been applied to turbulent flows (Hwang & Cossu 2010;McKeon & Sharma 2010).Particularly, the importance of choosing an appropriate eddy viscosity model was highlighted by Pickering et al. (2021).
In this study, we examine the intrinsic and extrinsic dynamics of STBLIs to reveal the nature of the low-frequency unsteadiness using GSA and resolvent analysis under a unified framework, in which low-frequency, large-scale motions are considered as coherent structures and turbulent fluctuations are closed by turbulence modelling.

Geometric configuration and flow conditions
The considered geometry is a compression corner with a ramp angle of α = 25°(see figure 1).A Cartesian coordinate system is constructed with the origin at the corner, the x direction along the flat plate, the y direction perpendicular to the flat plate and the z direction satisfying the right-hand rule.The flow conditions are taken from the experiments of Zheltovodov et al. (1990).The lowRe case has M ∞ = 2.95 and Re δ = 63 560 with δ = 2.27 mm measured 15.4δ upstream of the corner.The free-stream density is 0.314 kg m −3 , and the free-stream temperature is 108 K.The highRe case has M ∞ = 2.88 and Re δ = 132 840 with δ = 4.1 mm measured 8.04δ upstream of the corner.The free-stream density is 0.368 kg m −3 , and the free-stream temperature is 114.8K.The length of the flat plate is determined by matching the experimental boundary-layer thickness for the highRe and lowRe cases, respectively.

Governing equations
The flow is governed by the 3-D compressible RANS equations as where U is the vector of conservative variables and N is the nonlinear RANS operator.
Air is modelled as a calorically perfect gas with a specific heat ratio of 1.4 and a Prandtl number of 0.72.Sutherland's law is used to calculate the molecular viscosity.The Reynolds stresses are modelled using the Boussinesq assumption.The eddy viscosity is obtained using the Spalart-Allmaras (S-A) turbulence model (Spalart & Allmaras 1992) with the modification of Edwards & Chandra (1996).The reason why the S-A model is chosen is twofold.First, the length of the separation region predicted by the S-A model agrees reasonably well with existing experimental and LES data for the considered flow conditions as seen later.Second, the S-A model has been shown to be effective in predicting shock buffet, which is also a STBLI problem.The turbulent Prandtl number is set to 0.9.The flow-field variables are non-dimensionalized using the free-stream parameters and a characteristic length of L = 1 mm.

Flow solver
The 2-D and 3-D simulations are performed using an in-house finite-volume solver called PHAROS (Hao, Wang & Lee 2016;Hao & Wen 2020).The inviscid fluxes are calculated using the modified Steger-Warming scheme (MacCormack 2014), while the viscous fluxes are computed with a second-order central difference.An implicit line relaxation method (Wright, Candler & Bose 1998) is employed for pseudo time stepping, while a second-order implicit scheme is used for time-accurate simulations.The free-stream conditions are specified on the far-field boundary.A simple extrapolation is applied to the outflow boundary.The model surface is assumed to be no slip with a fixed wall temperature of 275.4 K. Computational grids with different levels of resolution are constructed to ensure grid independence (see Appendix A).Stability analysis requires a clean base flow with little numerical noise.For all cases, a reduction of nine orders of magnitude in the Euclidean norm of the density residual is achieved (see Appendix A).

Global stability analysis
It is assumed that vector U can be decomposed into a 2-D steady base flow (denoted by an overbar) and a small-amplitude 3-D unsteady perturbation (denoted by a prime) as U(x, y, z, t) = Ū(x, y) + U (x, y, z, t).
(3.2) Substituting (3.2) into (3.1) and neglecting the higher-order terms leads to the governing equations of the perturbation as where L is the linearized RANS operator.The expression of the linearized source term in the S-A model is given by Crouch et al. (2007).In other words, the turbulence model equation is also linearized without simplifications such as the frozen eddy-viscosity approach (cf.Carini et al. 2017).The perturbation U is further assumed to be in the following modal form: where Û is the 2-D eigenfunction, β is the spanwise wavenumber, ω r is the angular frequency and ω i is the growth rate.Substituting (3.4) into (3.3)leads to an eigenvalue problem, which is discretized in the same way as in the flow solver except that the inviscid fluxes are calculated using a central scheme in smooth regions detected by a shock sensor to reduce numerical dissipation.The boundary conditions are consistent with those for the base-flow simulation.Sponge layers are placed near the far-field and outflow boundaries to avoid any reflection of perturbations (Mani 2012).Note that the computational grids for the base-flow computation are directly used for the GSA with no data interpolation.The grid independence of the GSA is verified in Appendix A, where the effect of the size of the computational domain is also examined.
The discretized eigenvalue problem is solved using the implicitly restarted Arnoldi method implemented in ARPACK (Sorensen et al. 1996) for a given β in the shift-invert mode.The inversion is achieved using lower-upper decomposition implemented in SuperLU (Li et al. 1999).For each β, 50 eigenvalues are requested with a Krylov subspace of size 120.All the eigenvalues are converged with a residual less than 10 −9 .The flow is globally unstable if an eigenvalue with ω i > 0 can be found.An eigenmode is called stationary if ω r = 0 and oscillatory if ω r / = 0.The GSA solver has been applied to laminar interactions over various configurations (Hao et al. 2021(Hao et al. , 2022)).Here, it is extended to the RANS equations.

Resolvent analysis
To study the behaviour of the base flow as a noise amplifier, a small-amplitude forcing term f is added to (3.3) as where operator B constrains the forcing to a localized region and to the momentum components (Sartor et al. 2015;Bugeat et al. 2019).Both the forcing and response are assumed to be harmonic in time and in the spanwise direction as (3.7) Substituting (3.6) and (3.7) into (3.5)gives where R is the resolvent operator and I is the identity operator.The resolvent analysis seeks the forcing and response pair that maximizes the energy amplification or the optimal gain σ .The energies of the forcing and response are evaluated using the Euclidean norm and the Chu energy norm (Chu 1965), respectively.The optimization problem is converted to an eigenvalue problem (Bugeat et al. 2019), which is discretized in the same way as in the GSA and solved using the power iteration.

4.1.
Base flows Base-flow solutions are obtained using PHAROS for different Reynolds numbers and ramp angles.The streamwise velocity contours superimposed with the isoline of ū/u ∞ = 0.99 and the dividing streamline are displayed in figure 2 for the lowRe and highRe cases at α = 25°.Both cases exhibit large flow separation with the separation shock penetrating deeply into the incoming boundary layer.Figure 3 compares the 2-D RANS results with the experimental data of Zheltovodov et al. (1990) and the LES results of Grilli et al. (2012) in terms of the skin friction coefficient (C f = τ w /0.5ρ ∞ u 2 ∞ , where τ w is the wall friction) and wall pressure (normalized by the free-stream pressure p ∞ ) distributions.The skin friction drops suddenly to a negative value near the separation point on the flat plate and recovers near the reattachment point on the ramp.The pressure rises in the separation and reattachment processes, which are connected by a pressure plateau region.In general, the 2-D RANS simulation agrees well with the experiments and LES with respect to the separation location and the plateau pressure, except that the skin friction on the ramp is underestimated.

Global instability
In this section, we examine the global stability of the lowRe and highRe cases at α = 25°.Figure 4 presents the growth rate of the most unstable mode as a function of spanwise wavenumber.For the highRe case, the growth rate peaks at βL = 0.36, and the corresponding eigenvalue spectrum is shown in figure 5(a).The most unstable mode (referred to as the bubble mode) is stationary, whose streamwise and spanwise velocity perturbations (u and w ) are present in both the separation bubble and the reattached boundary layer (see figure 5c,d).The perturbation u has the same sign throughout the computational domain, while w has opposite signs in the upstream and downstream halves of the separation bubble.This mode greatly resembles the leading global mode in supersonic laminar interactions over a compression corner (Hao et al. 2023).
Underneath the most unstable mode, there is another stationary mode (see figure 5a).As β is decreased, these two modes move towards each other, leading to a modal coalescence that has also been reported in laminar interactions (Hao et al. 2021).As a result, a pair of At different ramp angles, the same modal features as for the 25°cases are obtained.As α is increased, the growth rate of the leading mode increases, and the corresponding spanwise wavenumber decreases (see figure 6).The flow becomes globally unstable as the ramp angle is beyond a critical value.Notably, the wavelength of the leading mode scales 971 A28-7 with the length of the separation region L sep , as shown in figure 6(b).Here, L sep is defined as the axial distance between the separation and reattachment points.
A 3-D unsteady RANS simulation is performed using PHAROS for the highRe case at α = 25°, which is initialized by duplicating the 2-D base flow in the spanwise direction.No initial disturbances are introduced, which means that global instability is expected to arise from the numerical round-off error.Only one spanwise wavelength of the most unstable mode is considered with 50 grid cells.The time step is set to 20 ns.Either increasing the number of spanwise grid cells or lowering the time step does not change the flow evolution (see Appendix A).The time history of the root mean square of the spanwise velocity (σ w ) on a wall-normal plane through x/L = 2.65 is shown in figure 7(a).After some initial transients, σ w experiences an exponential growth following the GSA prediction, which verifies the effectiveness of both methods.The contour of spanwise velocity on an x-y plane through z/L = 2.16 in the linear stage (tu ∞ /L = 6195) further confirms the emergence of the most unstable mode by comparing figures 7(b) and 5(d).Figure 7(c) shows the iso-surfaces of |w/u ∞ | = 1.5 × 10 −4 at tu ∞ /L = 6195 obtained from the 3-D RANS simulation.Also shown in the figure is the streamwise velocity contour at z/L = 0 to highlight the flow structure.A 3-D representation (see figure 7d) of the spanwise velocity perturbation is constructed using the eigenfunction of the bubble mode according to (3.4), which confirms that the exponential growth is governed by the bubble mode.Note that there is a phase difference in the spanwise direction between the RANS simulation and the GSA.Recall that the spanwise velocity perturbations have opposite signs in the upstream and downstream halves of the separation region, which stretch along the shear layer and the ramp surface, respectively.As seen later, such a perturbation field induces a pair of streamwise streaks in the reattached boundary layer after reaching a certain amplitude.
Following the exponential growth, the flow saturates into a new steady state.No supercritical Hopf bifurcation that leads to unsteadiness is found until the end of the simulation.Shown in figure 8(a-c), three wall-normal slices are extracted at x/L = 2.65, 20 and 40 in the saturated stage (tu ∞ /L = 29 734).The in-plane streamlines indicate a pair of counterrotating streamwise vortices along the ramp.The contour of the skin friction coefficient in figure 8(d) shows corrugated separation and reattachment lines in a zigzag pattern.Note that the contour is extended to two periods in the spanwise direction.As shown in figure 8(e), the presence of the streamwise vortices elevates the spanwise-averaged C f in comparison with the base-flow result, leading to better agreement with the experiments and the LES on the ramp.
In the 3-D simulation, the eddy viscosity is allowed to evolve.One may argue that the increasing eddy viscosity may have an unphysical stabilizing effect on the ending of the linear growth and any secondary instabilities after saturation.In this study, we fully rely on the predictive capability of the S-A model for 3-D large-scale flow structures.This is  consistent with previous unsteady RANS studies on shock buffet (Crouch et al. 2009;He & Timme 2021), where the eddy viscosity was unfrozen throughout the simulations.The effect of different RANS models on the flow bifurcation will be examined in a future study.

Response to upstream disturbances
We further examine the response of globally stable flows to external disturbances using resolvent analysis.Unless otherwise specified, the forcing is always localized approximately 10δ upstream of the corner to represent upstream disturbances, which extends from the wall to the far field.In fact, the optimal response is found to be insensitive to the forcing location (see Appendix C). Figure 9(a) shows the optimal-gain curves at ω r L/u ∞ = 10 −4 , 10 −3 , 10 −2 and 10 −1 as a function of spanwise wavenumber for the lowRe case at α = 23°.A low-pass feature similar to laminar interactions (Bugeat et al. 2022) is observed.Further lowering the frequency no longer changes the optimal-gain curve.At low frequencies, there are two local maxima in optimal gain at βL = 0 and 0.7.At high frequencies, the maximum optimal gain is found at βL = 1.2.The optimal gain and its frequency-premultiplied value are plotted in figure 9(b-d) as a function of angular frequency at βL = 0, 0.7 and 1.2.Only ω r L/u ∞ ≤ 0.2 is considered, which is at least one order of magnitude lower than the characteristic frequency in the incoming boundary layer and consistent with the assumption of separation of scales.However, this assumption holds less rigorously for the spanwise wavenumber.At the three wavenumbers, the weighted gain has a local maximum at ω r L/u ∞ = 6.5 × 10 −3 , 4.0 × 10 −4 and 7.5 × 10 −3 , respectively.At the three angular frequencies, the optimal response shown in the inset of figure 9(b-d) resembles the corresponding global mode predicted by the GSA, which indicates that the amplification is due to modal resonance.
Quantitatively, the 2-D optimal response at ω r L/u ∞ = 6.5 × 10 −3 is compared with the 2-D shock mode (ω r L/u ∞ = 0 and βL = 0) obtained from the GSA in figure 10(a   The global mode is supported by the separation bubble and is thus absent in the incoming boundary layer.The optimal response has a local peak at the forcing location, decreases monotonically along the flat plate and experiences an abrupt increase near the separation point.Inside the separation region, the Chu energy density of the optimal response agrees well with that of the global shock mode.The slight difference downstream of reattachment is likely due to the non-zero frequency of the optimal response and potential non-normalities.A similar comparison is made between the optimal response at ω r L/u ∞ = 4.0 × 10 −4 and βL = 0.7 and the bubble mode at the same spanwise wavenumber in figure 10(b).The almost overlapping Chu energy density distributions downstream of separation confirm that the local maxima of the optimal gain are essentially caused by modal resonance, while both convective-type non-normality (e.g. the Görtler instability) and component-type non-normality (e.g.transient growth of streamwise vortices) are insignificant (Chomaz 2005).
The optimal-forcing profiles corresponding to the optimal responses at ω r L/u ∞ = 6.5 × 10 −3 , 4.0 × 10 −4 and 7.5 × 10 −3 are plotted in figure 11.The 2-D optimal forcing works as a streamwise momentum pump in the incoming boundary layer, resulting in a periodic variation of the velocity profile.The 3-D forcings are characterized by streamwise vortices with the forcing energy contributed mostly by the spanwise momentum component.The spanwise component has its largest amplitude close to the wall and changes its sign near y/δ = 0.3.Interestingly, part of the forcing is present outside the boundary layer for both the 2-D and 3-D cases.Of particular interest is the 2-D optimal response as an excitation of the shock mode by external disturbances.Recall that the base flow is always stable to 2-D perturbations, which allows resolvent analysis at different Reynolds numbers and ramp angles.The resulting weighted gain is presented in figure 12.For the lowRe case, the angular frequency associated with the low-frequency peak decreases as α is increased.At α = 25°, increasing the Reynolds number further decreases the frequency.When non-dimensionalized using the length of the separation region, the frequencies with different interaction strengths collapse onto a universal value of fL sep /u ∞ = 0.015, which is consistent with the findings of Dussauge et al. (2006).
To reveal its nonlinear behaviour, a 2-D unsteady RANS simulation is performed, which is initialized using the 2-D base flow for the highRe case at α = 25°.The base flow is perturbed by the optimal forcing at fL sep /u ∞ = 0.015.The forcing amplitude is set to 0.01ρ ∞ u 2 ∞ /L.The time step is also 20 ns.The wall pressure is monitored at the base-flow separation point (x/L = −16.75),as shown in figure 13(a).The signal exhibits the well-known intermittent feature with the pressure oscillating between the free-stream pressure and the pressure behind the separation shock.The separation shock undergoes a back-and-forth motion (see the supplementary movie available at https://doi.org/10.1017/jfm.2023.687)with an intermittent region of approximately 0.6δ.The locations of the separation and reattachment points are shown in one cycle in figure 13(b).The bubble undergoes a breathing motion with the separation and reattachment points moving in opposite directions.

A28-12
On the low-frequency unsteadiness in STBLIs A more informative simulation is performed with the 2-D base flow perturbed by white noise at the same forcing location as the resolvent analysis for the highRe case at α = 25°.The random forcing is added to the right-hand side of the momentum equations with an amplitude of 0.1ρ ∞ u 2 ∞ /L and updated every 2500 steps to avoid violating the assumption of separation of scales.Consequently, the highest introduced frequency is 20 kHz (ω r L/u ∞ = 0.2).The simulation is run for 50 ms (tu ∞ /L = 30 973) with a time step of 20 ns. Figure 14(a) shows the wall pressure history at the base-flow separation point, which also has an intermittent behaviour.Spectral analysis is performed for the pressure signal between tu ∞ /L = 6195 and 30 973 with Welch's method (Welch 1967) to calculate the power spectral density (PSD) shown in figure 14(b).The signal is divided into three segments with 50% overlap.A Hann window is used for weighting the data.The frequency-premultiplied PSD peaks at fL sep /u ∞ = 0.015.This also suggests that the forcing profile is insignificant for modal resonance, which amplifies the low-frequency component of any upstream disturbances and results in a back-and-forth shock motion.
Additional RANS simulations are performed for the 3-D optimal responses at βL = 0.7 and 1.2 and other values of βL for the lowRe case at α = 23°.The forcing amplitude is set to 0.01ρ ∞ u 2 ∞ /L. Figure 15 presents the contours of C f at four instants with nearly equal intervals in one cycle at ω r L/u ∞ = 7.5 × 10 −3 and βL = 1.2.Note that the contours are extended to two periods in the spanwise direction.As expected, the saturated flow resembles that in § 4.2 featuring corrugated separation and reattachment lines and counterrotating streamwise vortices along the ramp.However, these vortices are not fixed in the spanwise direction but undergo a spanwise meandering motion (see the supplementary movie), which is reminiscent of the numerical results of Priebe et al. (2016) and Pasquariello et al. (2017).The contour of spanwise velocity on an x-y plane through z/L = 4.5 at tu ∞ /L = 1229 further confirms the excitation of the bubble mode identified by the GSA (see figure 16).

Conclusions
The low-frequency, large-scale motions in STBLIs over a compression corner are examined using GSA and resolvent analysis based on the assumption of separation of scales.
The GSA reveals that the 2-D base flow becomes unstable to 3-D perturbations as the ramp angle is beyond a critical value.Only a stationary unstable mode is identified, whose spanwise wavelength scales with the bubble size.At large spanwise wavenumbers, this mode is mostly confined to the separation bubble and the reattached boundary layer.At small spanwise wavenumbers, this mode manifests itself as a coupling between the separation shock and the separated shear layer.The GSA predictions are confirmed by a 3-D numerical simulation, where no supercritical Hopf bifurcation that leads to unsteadiness is found.The saturated flow exhibits a spanwise alternating expansion/contraction of the separation bubble and counterrotating streamwise vortices downstream of the corrugated reattachment line.For globally stable flows, the resolvent analysis identifies 2-D and 3-D local maxima in optimal gain due to modal resonance between external disturbances and the leading global mode, while the Görtler instability makes no evident contribution.A numerical simulation reveals that the 2-D optimal response is in the form of a large-scale, back-and-forth shock motion accompanied by a breathing bubble.The 3-D optimal responses all exhibit counterrotating streamwise vortices that meander in the spanwise direction with a range of frequencies and spanwise wavenumbers.This study suggests that the low-frequency shock motion in STBLIs is caused by the excitation of an intrinsic mode by extrinsic disturbances in line with the mathematical models of Plotkin (1975) and Touber & Sandham (2011), while the streamwise vortices are due to either global instability or modal resonance.is obtained, except that the k-ω model predicts an elevated skin friction on the ramp.The frozen eddy-viscosity approach (cf.Carini et al. 2017) is adopted in the GSA.In other words, the GSA solver for laminar flow (Hao et al. 2021) is used with the effective viscosity equal to the sum of the molecular and eddy viscosities.As a comparison, GSA is also performed on the base flow predicted by the S-A model using the frozen eddy-viscosity approach.For both cases, only a stationary unstable mode is captured.The growth rate curves as a function of βL are compared in figure 22.When the eddy viscosity is frozen, the modal coalescence phenomenon is absent, and the 2-D mode becomes globally unstable.This suggests that the 2-D unstable mode is likely an artefact of the frozen eddy-viscosity approach.Nonetheless, the obtained shapes of the leading mode at βL = 0 and 0.36 in figure 23 clearly exhibit the same features as those in figure 3, which indicates that the dynamics revealed in this study is insensitive to RANS modelling.

Figure 1 .
Figure1.Schematic of the geometric configuration and computational domain.

Figure 4 .Figure 5 .
Figure 4. Growth rates of the most unstable mode as a function of βL for different cases at α = 25°.The horizontal dashed line indicates zero growth rate.

Figure 6 .
Figure 6.(a) Growth rates and (b) spanwise wavenumbers and wavelengths of the most unstable bubble mode as a function of α for different cases.

Figure 7 Figure 8
Figure 7. (a) Time history of σ w at x/L = 2.65 in comparison with the growth rate of the most unstable bubble mode.(b) Spanwise velocity contour on an x-y plane through z/L = 2.16 at tu ∞ /L = 6195.(c) Iso-surfaces of |w/u ∞ | = 1.5 × 10 −4 and the streamwise velocity contour on an x-y plane through z/L = 0 at tu ∞ /L = 6195.(d) Iso-surfaces of |w /u ∞ | = 4 × 10 −7 reconstructed from the most unstable bubble mode at βL = 0.36 and the streamwise velocity contour obtained from the base-flow solution.Here Re δ = 132 840 and α = 25°.

Figure 9 .
Figure 9. (a) Optimal gains at different angular frequencies and (b-d) optimal gains and the frequency-premultiplied values at βL = 0, 0.7 and 1.2 and the corresponding optimal responses at Re δ = 63 560 and α = 23°.Local maxima are indicated by the vertical dashed lines and open circles.
) by plotting the distributions of the Chu energy density integrated in the wall-normal direction.The Chu energy densities are normalized by their respective maximum values.

Figure 10 .
Figure10.Distributions of the normalized Chu energy density integrated in the wall-normal direction for (a) the optimal response at ω r L/u ∞ = 6.5 × 10 −3 and βL = 0 and the shock mode at βL = 0 and (b) the optimal response at ω r L/u ∞ = 4.0 × 10 −4 and βL = 0.7 and the bubble mode at βL = 0.7 at Re δ = 63 560 and α = 23°.The vertical dashed lines mark the separation and reattachment points.

Figure 14
Figure 14.(a) Time history of wall pressure at the base-flow separation point and (b) the corresponding PSD for the highRe case at α = 25°.

Figure 19 Figure 20 .Figure 21 Figure 22 .Figure 23
Figure 19.(a) Schematic of two computational domains of different sizes and (b) eigenvalue spectra at βL = 0.36 obtained with different domains at Re δ = 132 840 and α = 25°.The dashed lines in (a) indicate the isoline of ū/u ∞ = 0.99 and the dividing streamline.The horizontal and vertical dashed lines in (b) indicate zero growth rate and zero angular frequency, respectively.