Low-frequency unsteadiness mechanisms in shock wave/turbulent boundary layer interactions over a backward-facing step

The low-frequency unsteady motions behind a backward-facing step (BFS) in a turbulent flow at $Ma=1.7$ and $Re_\infty=1.3718\times 10^5$ is investigated using a well-resolved large-eddy simulation (LES). The instantaneous flow field illustrates the unsteady phenomena of the shock wave/boundary layer interaction (SWBLI) system, including vortex shedding in the shear layer, the flapping motions of the shock and breathing of the separation bubble, streamwise streaks near the wall and arc-shaped vortices in the turbulent boundary layer downstream of the separation bubble. A spectral analysis reveals that the low-frequency behaviour of the system is related to the interaction between shock wave and separated shear layer, while the medium-frequency motions are associated with the shedding of shear layer vortices. Using a three-dimensional dynamic mode decomposition (DMD), we analyse the individual contributions of selected modes to the unsteadiness of the shock and streamwise-elongated vortices around the reattachment region. G\"ortler-like vortices, which are induced by the centrifugal forces originating from the strong curvature of the streamlines in the reattachment region, are strongly correlated with the low-frequency unsteadiness in the current BFS case. Our DMD analysis and the comparison with an identical but laminar case provide evidence that these unsteady G\"ortler-like vortices are affected by fluctuations in the incoming boundary layer. Compared to SWBLI in flat plate and ramp configurations, we observe a slightly higher non-dimensional frequency (based on the separation length) of the low-frequency mode.


Introduction
Shock wave/boundary layer interaction (SWBLI) has been an active research topic in the aerospace community over the past decades. This flow phenomenon is ubiquitous in high-speed aerodynamics, such as supersonic inlets, over-expanded nozzles, high-speed aerofoils (Green 1970;Dolling 2001). Shock-induced boundary-layer separation is a main contributor to flight drag of transonic aerofoils and pressure loss in engine inlets, which illustrates its relevance. Moreover, significant fluctuations of pressure and temperature are widely observed around the interaction regions. SWBLI can cause intense localized mechanical and thermal loads, which may eventually lead to the failure of material and structural integrity (Délery & Dussauge 2009;Gaitonde 2015). It is therefore crucial to take the effects of SWBLI into account in the process of aircraft design and maintenance, including material selection, assessment of fatigue life and thermal protection systems.
Canonical two-dimensional SWBLI configurations can be abstracted into three simplified cases: (1) incident (impinging-reflecting) shock, (2) compression ramp and (3) backward/forwardfacing step (BFS/FFS). Considerable progress has been achieved in understanding the unsteady phenomena and underlying mechanisms of SWBLI by means of advanced flow measurement techniques and well-resolved numerical simulations, particularly for the flat plate impinging shock and compression ramp configurations (Ganapathisubramani et al. 2007;Grilli et al. 2013;Pasquariello et al. 2017). These two cases share similar mean flow topology although the shocks are produced by different mechanisms, as shown in Figure 1(a) and (b). In the impinging/reflecting shock case, the incident shock induces a strong adverse pressure gradient on the boundary layer, which leads to the separation of the boundary layer. A separation shock is produced ahead of the separation point and a reattachment shock is generated around the reattachment location due to the compression of the boundary layer. For the ramp case, the strong flow compression caused by the ramp geometry induces a strong (separation) shock, which results in the separation of the incoming boundary layer. Subsequently, a reattachment shock is generated as the separated shear layer reattaches on the ramp downstream. In both cases, the SWBLI is accompanied by energetic unsteady motions at frequencies that are one or two orders lower than the boundary layer characteristic frequency ∞ / (Touber & Sandham 2009, 2011. Considerable research effort has been put into tracing the source of this low-frequency unsteadiness. In general, theories regarding the origin of this low-frequency motion of the separation shock are categorized as resulting from either upstream or downstream dynamics. The first group of theories associates the unsteady motions with upstream fluctuations within the incoming turbulent boundary layer. In an early work, Plotkin (1975) proposed a simple linear restoring model to explain the source of the shock wave oscillations, in which the shock is displaced by velocity fluctuations inside the upstream turbulent boundary and tends to return to its mean location through a restoring mechanism determined by the stability of the mean flow. The pressure measurement by Andreopoulos & Muck (1987) provided the first experimental evidence for a correlation of the shock wave unsteadiness with bursting events upstream the boundary layer in a compression ramp case at = 1.7. Unalmis & Dolling (1996) found low-frequency pressure fluctuations along the spanwise direction in the incoming boundary layer by measuring the pressure signal in the ramp case at = 5. Poggie & Smits (2001) performed measurements of wall pressure fluctuations and schlieren visualization in a backward-facing step/ramp configuration at = 2.9. They reported that also in this case the shock motion was correlated with upstream large-scale wave structures. Based on the cross-correlation analysis, they concluded that their experimental results are in good agreement with the linear restoring mechanisms proposed by Plotkin (1975). Beresh et al. (2002) used particle image velocimetry (PIV) and high-frequency response wall pressure transducers for a compression ramp interaction, and they found a clear correlation between streamwise velocity fluctuations in the lower part of the upstream boundary layer and low-frequency shock motions. In addition, they found no correlation between shock oscillations and the velocity fluctuations in the upper part of the upstream boundary layer, as well as the variation of the upstream boundary layer thickness, as reported by McClure (1992) in earlier work. Ganapathisubramani et al. (2007) also observed elongated superstructures with low-and high-speed streaks upstream the separation region in their stereoscopic PIV and planar laser scattering (PLS) measurements of a Mach 2 compression ramp interaction and they proposed these upstream large-scale structures are responsible for the low-frequency unsteadiness of the interaction region. Humble et al. (2009) further confirmed the presence of streamwise-elongated low-and high-speed streaks inside the upstream boundary layer using tomographic PIV for an incident shock interaction at = 2.1. Their results show that this reorganization of the upstream boundary layer in both streamwise and spanwise directions conforms to the overall streamwise translation and spanwise rippling of the interaction region. However, Touber & Sandham (2011) argued that the low-frequency interaction motions do not necessarily require a forcing source from upstream or downstream and are more like an intrinsic response to the broadband frequency spectrum of the upstream turbulent fluctuations. Porter & Poggie (2019) consider that this response is a selective response of the separation region to certain large-scale perturbations in the lower half part of the upstream boundary layer based on their high-fidelity simulation.
The second group of theories attributes the low-frequency dynamics to mechanisms intrinsic to the interaction system itself, that is, with an origin downstream of the separation line. Already early experimental studies suggested that the low-frequency motion of the separation shock is linked to the expansion and contraction of the separation bubble (Erengil & Dolling 1991;Thomas et al. 1994). For the impinging shock induced interaction, Dupont et al. (2006) found a clear statistical link between low-frequency oscillation of the separation shock and the downstream interaction region by analysing experimental pressure signals. Furthermore, they also reported a quasi-linear relation between the separation shock and the reattachment shock motions. By DNS of a Mach 2.25 impinging shock case, Pirozzoli & Grasso (2006) established a resonance theory, in which acoustic waves are produced by the interaction between coherent structures in the bubble and the incident shock. The upstream propagation of these acoustic waves is responsible for the low-frequency oscillations of the SWBLI system. Touber & Sandham (2009) performed a global linear stability analysis of the mean flow field from their LES and detected an unstable global mode inside the separation bubble, which provides a possible driving mechanism for the low-frequency unsteadiness by displacing the separation and reattachment points. Piponniau et al. (2009) proposed a simple physical model that relates the low-frequency oscillations to the breathing motions of the separation bubble, in which the collapse of the separation bubble is caused by a continuous entrainment of mass flux, while the dilation corresponds to a radical expulsion of the mass injection in the bubble. A similar model was suggested by Wu & Martín (2008) based on DNS of a compression ramp configuration. They consider that a feedback loop, involving the separation bubble, the detached shear layer and the shock system, is the underlying mechanism for low-frequency shock motions. The DMD analysis of Grilli et al. (2012) provided further evidence that mixing across the separated shear layer leading to a contraction and expansion of the separation bubble is the dominant mechanism for the low-frequency unsteadiness. Numerical work of Grilli et al. (2013) and Priebe et al. (2016) identified streamwise-elongated Görtler vortices originating around the reattachment location for compression ramp configurations. For an impinging shock configuration, Pasquariello et al. (2017) reported very similar observations of low-frequency DMD modes characterised by streamwise-elongated regions of low and high momentum that are induced through Görtler-like vortices. As the separation-bubble dynamics is clearly coupled to these vortices, Görtler-like vortices might act as a source for continuous (coherent) forcing of the separation-shock-system dynamics.
In an attempt to resolve this discrepancy, Souverein et al. (2010) proposed that actually both upstream and downstream mechanisms contribute to the SWBLI dynamics with case dependent intensity. Which type of mechanism is more dominant in producing the low-frequency dynamics depends on the shock strength and possibly the Reynolds number. In weak interactions, the low-frequency unsteady motions can be mainly associated with upstream effects, while the unsteadiness of the strong interactions are more likely driven by the dynamics of the downstream separation bubble and reattachment shock (Clemens & Narayanaswamy 2014). Also Priebe et al. (2016) implied that upstream disturbances contribute to the low-frequency behaviour although they consider that the downstream Görtler instability is the dominant one. Bonne et al.
(2019) indicated that the low-frequency oscillations involve both the amplification of upstream disturbances by the separated shear layer and a feedback excitation from the shock foot and backward travelling density waves.
As discussed above, SWBLI in the impinging shock and compression ramp configuration share similar unsteady behaviour and physical mechanisms (Clemens & Narayanaswamy 2014;Smits & Dussauge 2006). In contrast to these well-analysed canonical cases, supersonic flow over a BFS has a distinctly different flow topology, as shown in Figure 1(c). The incoming turbulent flow undergoes first a centred Prandtl-Meyer expansion (PME) with the separation location fixed at the step convex corner. The free shear layer then develops towards the downstream wall on which the flow reattaches. Compression waves are generated around the reattachment location, which coalesce into a reattachment shock (Loth et al. 1992;Sriram & Chakraborty 2011). In this configuration, the upstream limit of the separation bubble is stationary and only the downstream reattachment shock is present. The dynamics of the recirculation and shock region is reported to be unsteady as in other conventional cases (Bolgar et al. 2018). In an early experimental study, by examining the variation of skin friction, Ginoux (1971) observed the systematic development of counter-rotating streamwise vortices around the reattachment, occurring in laminar, transitional and turbulent flows alike. The wavelength of these vortices is equal to two or three times the boundary layer thickness for a wide range of Mach numbers. These Görtler-like vortices were also reported in the experimental visualization via nano-tracer-based planar laser scattering (NPLS) (Zhu et al. 2015). In addition, small unsteady shedding vortices along the shear layer were identified by Chen et al. (2012) using the same visualization techniques. However, the Kelvin-Helmholtz (K-H) vortices typical in laminar and transitional cases were not observable in the turbulent shear layer (Zhi et al. 2014). The observed coherent vortical structures cover a wide range of length and frequency scales, involving the vortex shedding close to the step, longitudinal vortices and hairpin vortices downstream of the shear layer (Soni et al. 2017). The unsteady characteristics can be quantified by the dimensionless Strouhal number = / ∞ based on the reattachment length and free stream velocity. By means of particle image velocimetry and dynamic pressure measurements, Bolgar et al. (2018) inferred that for a flow at = 2.0 the higher frequency content ( = 0.05 − 0.2) is related to the shock motions, while the dominant low-frequency parts ( ≈ 0.03) are associated with the separation bubble. More efforts are required to scrutinize the frequency characteristics of BFS SWBLI and to analyse whether the low-frequency unsteadiness of supersonic BFS flows has a similar origin as that in the impinging shock and ramp SWBLI cases. In our previous work (Hu et al. 2019(Hu et al. , 2020, we examined the unsteady SWBLI over a BFS in a laminar inflow regime. The preceding discussion motivates to investigate to what extent the laminar and turbulent cases share similar unsteady features and physical mechanisms. In this paper, we analyse new large-eddy simulation (LES) results for a fully turbulent BFS flow at = 1.7 with special attention to the low-frequency dynamics. The organization of the paper is as follows. Details of the numerical methods used and the setup of the flow configuration are given in §2. Then the flow topology of the mean and instantaneous flow is discussed in §3. The characteristic frequencies of the significant unsteady motions are analysed using spectral analysis. Finally, dominant modes in the SWBLI are extracted via a three-dimensional dynamic mode decomposition (DMD). By comparing with previous works, a physical mechanism of the low-frequency unsteadiness source is proposed ( §4). The conclusions with a summary of the main results are presented in §5.

Governing equations
The physical problem is governed by the unsteady three-dimensional compressible Navier-Stokes equations with appropriate boundary and initial conditions, and the constitutive relations for an ideal gas. We solve the conservation equations for mass, momentum and total energy 2) where is the density, the pressure and the velocity vector. The total energy is defined as (2.6) The fluid is assumed to behave as a perfect gas with a specific heat ratio = 1.4 and a specific gas constant = 287.05 J(kg · K) −1 , following the ideal-gas equation of state = . (2.7) The dynamic viscosity and thermal conductivity are a function of the static temperature and are modelled according to Sutherland's law and the assumption of a constant Prandtl number (2.8) The values adopted for the computations are: = 18.21 × 10 −6 Pa · s, = 293.15 K, = 110.4 K and = 0.72.

Flow configuration
The current computational case is an open BFS (i.e., no upper wall) with a supersonic turbulent boundary layer inflow, a schematic of which is shown in Figure 2. The origin of the Cartesian coordinate system is placed at the step corner. The turbulent inflow is characterised by the freestream Mach number ∞ = 1.7 and the Reynolds number 0 = 13718 based on the inlet boundary layer thickness 0 (99% ∞ ) and free-stream viscosity. The main flow parameters are summarized in Table 1   exclude potential uncertain effects from the numerical inlet boundary conditions on the flow in the region of interest. The height of the step ℎ = 3 0 is three times larger than the inlet boundary layer thickness. In addition to this (fully) turbulent BFS flow, we also present selected results for a fully laminar inflow case with the same free stream flow parameters and geometry (Hu et al. 2019) for comparison. Note that this laminar inflow case is referred to as the laminar case for the simplicity of the discussion although flow transition to turbulent occurs shortly downstream of the step.

Numerical method
The LES method of Hickel et al. (2014) is used to solve the governing equations. The sub-grid scale model is fully merged into a non-linear finite-volume scheme provided by the adaptive local description method (ALDM) (Hickel et al. 2006(Hickel et al. , 2014. ALDM is based on a solution-adaptive reconstruction operator and a numerical flux function that incorporates the essential elements of LES, filtering and deconvolution. The optimization procedure starts from a nonlinearly stable numerical scheme and towards a final ALDM scheme which acts as an accurate sub-grid scale model. The viscous flux is discretized by a second-order central difference scheme and an explicit third-order total variation diminishing (TVD) Runge-Kutta scheme (Gottlieb & Shu 1998) is used for time marching. This method has been successfully applied to various supersonic flow cases, including shock wave/boundary layer interactions (SWBLI) on a flat plate (Pasquariello et al. 2017) and compression ramp (Grilli et al. 2012(Grilli et al. , 2013, and transition between regular and irregular shock patterns in SWBLI (Matheis & Hickel 2015), as well as in our previous work on SWBLI in laminar and transitional BFS flows (Hu et al. 2019(Hu et al. , 2020. More details about the numerical method can be found in the literature (Hickel et al. 2006(Hickel et al. , 2014. The numerical grids are generated using a Cartesian grid structure with block-based local refinement, as displayed in Figure 3. In addition, hyperbolic grid stretching was used in the wall-normal direction downstream of the step. The mesh is sufficiently refined near all walls with + < 0.9 to ensure a well-resolved wall shear stress. The grid spacing becomes coarser with increasing wall distance but the expansion ratio between the adjacent blocks is not larger than two. The distribution of mesh cells are uniform in the spanwise direction. Using this discretization strategy, the computation domain has around 36 × 10 6 grid points and a spatial resolution of the flow field with Δ + max × Δ + max × Δ + max = 36 × 0.9 × 18 in wall units for the entire domain (Δ + max = 0.9 on the step wall). The temporal resolution, that is the time step, is approximately Δ ∞ / 0 = 7.6 × 10 −4 , corresponding to a Courant-Friedrichs-Lewy condition CFL 0.5.
The step and wall are modelled as no-slip adiabatic surfaces. All the flow variables are extrapolated at the outlet of the domain. On the top of the domain, non-reflecting boundary conditions based on Riemann invariants are used. Periodic boundary conditions are imposed in the spanwise direction. Inlet turbulent boundary conditions require a special approach since a very large domain for the natural development of turbulence is undesirable in view of computational resources and time. We use a synthetic turbulence generation method based on digital filter technique (Klein et al. 2003) to produce the appropriate turbulent inflow. This method can reproduce both first-and second-order statistical moments and spectra, without introducing lowfrequency content which may modulate the low-frequency dynamics downstream. The reference data used are from Petrache et al. (2011) to specify realistic integral length scales and mean boundary layer profiles. According to previous studies (Grilli et al. 2013;Wang et al. 2015), a transient length of around 10 0 is sufficient for turbulence to develop in the supersonic boundary layer under these conditions. Nevertheless, we place the inflow plane 40 0 upstream of the step.
In order to examine the grid and domain size independence, the van Driest transformed mean velocity profile and Reynolds stresses in Morkovin scaling are provided at / 0 = −5.0 in Figure 4. The computed flow field reached a fully developed statistically steady state after an initial transient period of ∞ / 0 = 800. The samples then were collected every ∞ / 0 = 0.5 over an interval of another ∞ / 0 = 600, yielding an ensemble size of 1200. For comparison, the figure also includes the theoretical law of the wall and incompressible DNS data of Schlatter & Örlü (2010) at = 360 and = 1000. The present mean velocity profile is consistent with both the logarithmic law of the wall (with the constants = 0.41 and = 5.2) and the DNS data. The Reynolds stresses from the current LES are also in a good agreement with the reference data. Since the current LES data is for a compressible boundary layer that has a higher momentum thickness Reynolds number = 2000 and friction Reynolds number = 400, the velocity profile has a slight larger plateau value and streamwise Reynolds stress profile features with a higher peak value in the buffer layer (Marxen & Zaki 2019).  The mean reattachment length (equal to = = 8.9 0 ≈ 3.0ℎ) is defined by the location of zero mean skin friction, = 0, in Figure 6(a). The value of increases upstream of the step due to the flow acceleration induced by the expansion near the separation point ( = 0). Behind the step, there is a 'dead-air' zone where the recirculating velocity is extremely low. Thus, uniform ≈ 0 are observed in the first 30% of the separation bubble (0.0 / 0 2.8). The separated flow then rapidly reaches its strong level at ≈ 2.1ℎ ≈ 6.2 0 , which is very close to the value ( ≈ 2ℎ ≈ 6.4 0 ) reported by Chakravarthy et al. (2018). As the free shear layer reattaches on the downstream wall ( / 0 = 8.9), the turbulent boundary layer recovers and returns to a typical turbulent level ( = 0.0027). The reattachment length ≈ 3.0ℎ based on the step height is in a good agreement with the previous experimental work by Bolgar et al. (2018) and the numerical study by Chakravarthy et al. (2018), who reported values of = 3.2ℎ and = 3.0ℎ, respectively. Compared with the laminar case (blue dotted lines), the mean skin friction further confirms the shorter separation length in the turbulent case. The turbulent case has a much higher upstream of the step. The laminar case reaches, however, a similar level downstream of the separation region, because laminar-to-turbulent transition is triggered within the separated shear layer. Figure 6(b) shows the streamwise variation of the wall pressure. As we can see, upstream the step, the wall pressure remains at almost the same level. The pressure drops drastically to around 42% ∞ in the first half of the separation bubble due to the expansion and the less energetic recirculating flow. The wall pressure then continues decreasing slowly to its global minimum at / 0 = 4.6, corresponding to the relatively strong reversed flow in terms of in Figure 6(a). As the boundary layer reattaches on the wall and undergoes compression, the wall pressure quickly returns to the initial level. In the experimental work of Hartfield et al. (1993), they reported that the measured pressure deceases from 34.8 kPa to around 14.2 kPa (≈ 41% ∞ ) upstream of the seperation bubble and returns to the free stream level downstream the interaction region, which are in a good agreement with the current results. In the laminar regime, the expansion fan is not as strong as for the turbulent case. Similarly, the intensity of the reattachment shock is weaker in the laminar case corresponding to a slower wall-pressure rise downstream. We see the expected small-scale coherent structures in the incoming turbulent boundary layer. Since the separated shear layer is inviscidly unstable, it rolls up and then larger and stronger vortical structures are generated over the bubble region. As the shear layer evolves downstream, the upstream small turbulent structures develop into larger coherent structures due to the shear layer instability, indicated by the arc-shaped vortices in the outer region of the boundary layer downstream the bubble. These coherent vortical structures propagate above the reversed flow from the separation to the reattachment location, and they also exist within the turbulent boundary layer downstream of the bubble.

Instantaneous flow organisation
For comparison, the instantaneous vortical structures of the laminar case are provided in Figure 7(b). The typical K-H vortex structure present in the laminar case is not observed in the current turbulent regime where the quasi two-dimensional vortices are probably distorted by the highly three-dimensional turbulence. In the middle of the shear layer, large coherent Λ-shaped vortices are formed and transformed into arc-shaped vortices downstream in the laminar case as a result of vortex stretching and tilting, whereas only arc-shaped vortices are present downstream in the turbulent case. From the numerical schlieren image shown on the / 0 = −4 slice, the shock intensity in the laminar case is weaker than that of the turbulent one, which is consistent with the evolution of the streamwise wall pressure in Figure 6(b).

Unsteady characteristics
The flow field over the BFS is highly unsteady, with vortices of various spatial scales observed in the visualization of Figure 7. To characterize the regions of most prominent unsteadiness, the variance of the wall-normal velocity is provided in Figure 8. As we can see, the most active region can be found along the separated shear layer (between the isoline of = 0 and boundary layer edge), especially in the proximity of the reattachment location with a maximum of approximately 0.18 ∞ occurring at / 0 = 7.2, / 0 = −2.2. These major fluctuations caused by the recompression have also been reported in previous experimental work ( Bolgar et al. 2018). Additionally, relatively weak fluctuations are found along the reattachment shock, reflecting its unsteady position. For the other normal Reynolds stress components and , high levels of fluctuations are similarly observed around the reattachment point. We see that the separated shear layer and shock wave system is highly unsteady over the BFS with similar fluctuation intensities as in other canonical SWBLI geometries (Touber & Sandham 2011;Pasquariello et al. 2017).
Our attention then is put on the zones of the shear layer, reattachment location and shock wave to scrutinize the dynamic motions by examining a number of snapshots of the instantaneous flow field. First of all, we take a closer look at the shear layer. Figure 9 displays the contours of the streamwise velocity and streamlines at two arbitrarily selected instants. There are positive and negative streamwise velocity fluctuations alternating along the shear layer, which is the expected footprint of the shear layer instability behind the step. The convective Mach number , defined as is ≈ 0.93 at / 0 = 4.5, where 1 and 2 are the maximum streamwise velocity at the high-speed and low-speed sides of the mixing layer, and 1 , 2 are the speed of sound at the corresponding locations. As indicated by Sandham & Reynolds (1991), the compressible shear layer also exhibit three-dimensional features at this convective Mach number, which explains the emergence of the three-dimensional waves in the shear layer behind the step shown in Figure 7. As the free shear flow evolves downstream, large coherent vortices are caused by the vortex  Figure 10(b) provides the spanwise wavenumber weighted power spectral density (PSD) of the streamwise wall velocity at two stations. As we can see, the wavenumber of the upstream structures ( / 0 = −5.0) is ≈ 2.0, corresponding to a spanwise wavelength ≈ 0.5 0 . The shear stress is relatively uniform at a low level downstream the step (0 < / 0 < 5.0) due to the less energetic flow in this region. Shortly upstream of the mean reattachment location (5 < / 0 < 8.9), there is significant reverse flow, cf. Figure 6(a), and indicates an increased spanwise length of the coherent structures. After reattachment, streamwise-oriented features are observed in the skin friction maps that indicate large scale streaks with a spanwise alternation of high and low velocity. For example at / 0 = 10.0, the dominant spanwise wavenumber of the streamwise wall velocity is ≈ 0.35 ( ≈ 2.9 0 ), as shown in Figure 10(b). Further downstream, the intensity of becomes more homogeneous again. Similar phenomena have been reported in previous experiments of BFS with a wide range of Mach number (Ginoux 1971). The up-wash and down-wash effects of the Görtler-like vortices are believed to induce the alternating low and high skin friction in the spanwise direction around the reattachment, as will be discussed in the following sections. The characteristic wavelength of these streaks is between = 2.0 0 and 3.3 0 , which is consistent with previous experimental and numerical observations, reporting that the wavelength of these vortices is between two and three times the boundary layer thickness (Ginoux 1971;Priebe & Martín 2012;Grilli et al. 2013;Pasquariello et al. 2017).
In addition to the these relatively local phenomena, a large-scale unsteady motion is identified in the interaction system, as shown by the instantaneous velocity fields at two instants in Figure 11. These two instants represent different states of the separation bubble, i.e., expansion and shrinking. At ∞ / 0 = 954.5, the length of separation bubble is around / 0 = 7.5, while the flow reattaches further downstream at about / 0 = 9.0 when expanding at ∞ / 0 = 1080. In addition, the position of the shock (marked as white isolines of |∇ | 0 / ∞ = 0.4) moves, most notably in the shock foot region. At ∞ / 0 = 954.5, the shock foot locates somewhere between / 0 = 7.5 ∼ 10.0 and the shock angle is = 22.2 • . At ∞ / 0 = 1080, shock foot is between / 0 = 5.0 ∼ 7.5 and shock angle reduces to = 16.8 • . It is clear from this comparison that the recirculation area and shock location vary in time.
For the laminar case, we also observe vortex shedding along the shear layer and the flapping   (Hu et al. 2019). However, there are notable differences in the near wall dynamics, as can be seen when comparing the instantaneous skin friction contours and spanwise wavenumber weighted power spectral density in Figure 10 (turbulent case) to 12 (laminar case). The distribution of the laminar case skin friction is obviously spanwise uniform upstream the step. As the separated shear layer undergoes laminar-to-turbulent transition, the skin friction contours develop weak two-dimensional features around the reattachment location and further downstream. The dominant spanwise wavenumber of the upstream ( / 0 = −5.0) and downstream ( / 0 = 10.0) structures is close to each other with ≈ 0.8. The large low-and high-speed due to the small spanwise wavelength ≈ 1.2 0 ( ≈ 2.9 0 around the reattachment in the turbulent case). This difference suggests that there is probably no evidence of the counter-rotating Görtler vortices in the laminar case.

Spectral Analysis
An overview of frequency characteristics for the shock wave and separated boundary layer system is provided by the frequency weighted power spectral density of the wall pressure at selected streamwise locations in Figure 13. The sampling interval is ∞ / 0 = 950 ∼ 1350 with a sample frequency 0 / ∞ = 4, excluding the initial transient stage of the simulation. Welch's method with Hanning window was applied to compute the PSD using eight segments with 50% overlap (the same for the following PSD calculations). Upstream of the step ( / 0 = −3.0), the spectrum shows a broadband bump centred around = 0 / ∞ = 0.8, which is close to the characteristic frequency ( ∞ / ) of the upstream turbulent boundary layer (Dolling 2001). The low-frequency contents are relatively small upstream of the step, which demonstrates that the digital filter technique does not introduce significant spurious low-frequency features into the boundary layer. Downstream the step, we observe broadband low-frequency content between = 0.01 ∼ 0.8 ( ℎ = ℎ/ ∞ = 0.03 ∼ 2.4), in addition to the typical signature of boundary layer turbulence at the higher frequencies. Two significant low frequencies can be identified along the streamwise distance. The lower one is around = 0.013 (lower blue dashed line in the graph), which is most significant in a short distance behind the step ( / 0 3.0). It appears that this low frequency is not the dominant one further downstream the separation bubble and an intermediate frequency at = 0.1 ∼ 0.3 (upper region separating by green dashed lines) begins to take the lead up to / 0 = 20.0. In the traditional ramp and impinging shock cases (Ganapathisubramani et al. 2007;Agostini et al. 2012;Pasquariello et al. 2017), the medium-frequency shear layer oscillations arise after the separation and the downstream propagation of this dynamics affects the reflected shock featuring intermediate frequencies, while the interaction between separation shock and boundary layer exhibits the low-frequency behaviour. Based on these conclusions, the medium frequency motions of the present BFS case are probably related to the shear layer instability, the downstream convection of which produces a relatively significant medium-frequency unsteadiness around the reattachment location ( / 0 = 9.25). For the low-frequency contents of the BFS case, they are likely connected to the interactions of the reattachment shock and the separation bubble, the feedback of which leads to the low-frequency peak immediately downstream of the step ( / 0 = 1.0).
To further confirm this conjecture, several aerodynamic parameters are extracted from the medium low Figure 13: Frequency weighted power spectral density of the wall pressure with the streamwise distance.
current results. For the medium-frequency behaviour, the temporal variation of the streamwise velocity within the separated shear layer and the spanwise-averaged reattachment position are plotted in Figure 14. These data are extracted with the same sampling frequency as the aforementioned pressure signal. The location of the spanwise-averaged reattachment point is obtained as follows: the isolines of the streamwise velocity = 0 are collected at each time step; and in each spanwise plane the most downstream position meeting this condition ( = 0) is determined as the instantaneous value of . An unsteady motion at a frequency around = 0.2 ( ℎ = 0.6) appears energetically dominant for both shear layer velocity and reattachment location, which is more clear in the spectra of Figure 14. This medium frequency is the characteristic frequency of the shedding vortices within the shear layer. These vortices are shedding downstream as the shear layer and pass through the reattachment downstream the bubble, which explains that a similar frequency is observed in the spectrum of the reattachment location. There are also less energetic peaks at lower frequencies around = 0.03, which will be discussed in the next paragraph. When taking a closer look on a short interval in Figure 15, the velocity signal of the shear layer is more periodic and regular. In contrast, the curve for the reattachment point follows a more sawtooth-like trajectory, along which its value undergoes a sharp drop when the reattachment point moves upstream, while it experiences a less rapid relaxation as the reattachment location shifts downstream, for instance around ∞ / 0 = 1160. The sawtooth-like behaviour was also reported for incident shock and ramp cases (Priebe & Martín 2012;Pasquariello et al. 2017), and is attributed to the passage of shedding vortices formed in the shear layer near the reattachment.
With regard to the global dynamics, the temporal variation of the spanwise-averaged reattachment shock angle and separation bubble volume are shown in Figure 16. The bubble volume per unit spanwise length is defined as the area between the isoline of = 0 and the bottom wall. The shock angle is determined based on the pressure gradient outside the boundary layer by fitting the isolines of |∇ | 0 / ∞ = 0.24. We obtain two values by intersecting the isolines of |∇ | 0 / ∞ = 0.24 at / 0 = 0.5 and then take the average of these two values as the first streamwise coordinate of the shock position. A second point of the shock position is obtained by repeating the same operation at / 0 = 5.0. A straight line is fitted based on these two points and the angle between the fitting line and the −direction is considered as the shock angle. Both curves of the separation bubble size and shock angle are irregular and aperiodic in time, which suggests that the unsteady motion involves a range of time scales, cf. Refs (Dussauge et al. 2006;Priebe et al. 2016). For the signal of the separation bubble volume, shown in Figure 16(a), there is a significant low-frequency peak at = 0.023 in the spectrum. It indicates that the bubble expands and shrinks with a frequency whose value is about two-order lower than the frequency of the typical turbulence. The spectrum of the shock angle also displays a peak at = 0.023,  see Figure 16(b), which is much more pronounced than the peak observed for the reattachment location at the same frequency in Figure 14(b). In addition, there is a second frequency peak around = 0.13, which corresponds to the dominant frequency in the spectrum of reattachment location. Since the shock is formed by the compression waves originating at reattachment, spectra for the shock and reattachment locations include peaks at common frequencies.
The statistical connection between the low-frequency signals can be quantified through coherence and phase . The spectral coherence between two temporal signals ( ) and ( ) is defined as where is the power spectral density of ( ) and ( ) represents the cross-power spectral density between signals ( ) and ( ). The phase is determined by For a specific frequency, if 0 < < 1, it means that there is noise in the datasets or the relation between these two signals is not linear. When equals to 1, it indicates that the signals ( ) and ( ) are linearly related, and = 0 signifies that they are completely unrelated. The coherence and phase between the separation bubble volume and shock location of the spanwise-averaged snapshots are shown in Figure 17. The definition of the separation bubble volume is the same as before. The shock location is the coordinate of the intersection between -axis and the fitted straight shock line (defined when calculating the shock angle in Figure 16). A high value of coherence ( = 0.42) is observed at the frequency = 0.028, which manifests that the separation bubble and reattachment shock are nonlinearly related to each other around the shown dominant low frequency in the spectrum of Figure 16. Moreover, these two signals are approximately in phase, as can be seen from the low level of . The above observations provide evidence that the unsteady low-frequency behaviour is related to the breathing of the separation bubble and the flapping motion of the shock, while the medium-frequency motions are associated with the shedding vortices of the shear layer. Thus a decoupling of the frequency scales is required to further trace the sustained source of the intrinsic unsteadiness of the interaction, which is the objective of §3.5.
Similar low-and medium-frequency are also observed for the laminar cases. Figure 18 plots the corresponding frequency weighted power spectral density of streamwise velocity around the mean reattachment location and the spanwise-averaged separation bubble size. To compare the laminar and turbulent cases, the frequency is rescaled by the reattachment length as = / ∞ . For the signal of streamwise velocity in Figure 18(a), the results show a broadband low-frequency spectrum for both the laminar and turbulence cases. However, a local spectrum peak is clearly observed at = 0.15 ( = 0.017) in the turbulent case and at = 0.20 ( = 0.018) for the laminar case. Since there are distinct shedding vortices in the shear layer for both flow regimes, the relevant prevailing medium frequencies are close to each other. For the bubble size in   Figure 18(b), the dominant frequency of the separation bubble in the laminar case is = 0.33, while the corresponding value is lower ( = 0.22) in the turbulent case. These differences suggest that there are probably other flow dynamics involved, which leads to a lower frequency of the unsteady motions in the turbulent case. As previously discussed, Görtler-like vortices are likely to be associated with the low-frequency unsteadiness of SWBLI (Priebe et al. 2016). Therefore, we infer that the streamwise streaks in the turbulent regime may play a role in the transformation of the dominant low frequency.

DMD analysis of the three-dimensional flow field
To better separate the different dynamics from the coupled broadband frequency spectrum, a frequency-orthogonal modal decomposition of the three-dimensional flow field is conducted based on dynamic mode decomposition (DMD). Schmid (2010) first proposed this method to identify the most important dynamic information contained in equal-interval temporal snapshots from an unsteady flow field. Briefly summarized, a set of (reduced) modes will be extracted from the original dynamic system, each of which is associated with a single frequency behaviour and the combination of which approximate the complete unsteady system. Compared with proper orthogonal decomposition (POD), which is usually used for obtaining a low-dimensional reconstruction of the dynamic system, DMD focuses on the relevant flow dynamics while decoupling in frequency. This technique has been widely applied for various unsteady flow problems, such as the transition mechanism from laminar to turbulent flow (Hu et al. 2019), unsteadiness of the SWBLI (Pasquariello et al. 2017) and the identification of coherent vortex structures (Wang et al. 2020).
Following the DMD methodology, the original dynamic system can be represented by ( 3.4) where can be considered as the amplitude of the -th DMD mode and the Vandermonde matrix and signifies the temporal evolution of the dynamic modes. The eigenvalues are usually further converted into a more familiar complex stability plane through the logarithmic mapping = ln( )/Δ (Leroux & Cordier 2016). The dynamic information about the growth rate and angular frequency of a specific DMD mode are then computed by To facilitate a physical interpretation, we also reconstructed the real-valued flow field of the individual modes by superimposing the fluctuations from each mode onto the mean flow , formulated as where and are the amplitude and optional amplification factor of the corresponding mode . The reconstructed flow field at different phase angles represents the temporal evolution of the dynamic system, In this way, the imaginary part of the reconstruction at a phase angle = 0 is equivalent to the real part at = , and vice versa.
In the above analysis, we identified two types of unsteady behaviour at different frequencies. However, part of the signals were extracted from the spanwise-averaged field, like reattachment location, bubble size and shock angle; thus spanwise unsteady features may be missing from the two-dimensional flow field and a three-dimensional DMD analysis is required. Considering the large size of the data ensemble, a spatial subdomain (−5.0 / 0 20.0 and −3.0 / 0 5.0, covering the most interesting region) is extracted with a downsampled spatial resolution in all directions. The present DMD analysis of the three-dimensional subdomain is carried out based on 1200 equal-interval snapshots with the same temporal range of the previous signals and a smaller sampling frequency 0 / ∞ = 2 as the frequencies above the characteristic frequency of the turbulent integral scale ∞ / 0 are not of our current interest. The resulting frequency resolution is 1.67 × 10 −3 1. Figure 19(a) shows the calculated eigenvalue spectrum provided by the standard DMD. The obtained modes appear as complex conjugate pairs and most of them are well distributed along the unit circle | | = 1 except a few decaying modes within the circle, which means the resulting modes are saturated (Rowley et al. 2009). The magnitudes of the normalized amplitudes (| | = | |/| | max ) of the corresponding DMD modes are shown in Figure 19(b) for the positive frequencies and gray shaded by the growth rate . Here, the strongly decaying modes (| | 0.95) have been removed, as they do not contribute to the long-time flow evolution. The darker the vertical lines are, the less decayed the modes are. The convergence of the DMD  Table 2: Information of the selected modes results was verified by repeating the DMD using 400 snapshots less, which confirmed that the current DMD results are well-converged with respect to the number of the snapshots. From the frequency-magnitude spectrum, we identified three interesting frequencies, a lower one (marked as A) with < 0.06, a medium one (marked as B) with 0.06 0.2, and a higher one (marked as C) with > 0.2. Based on the growth rate and magnitudes of the modes, three modes are selected from the frequency spectrum, one representative for each of the frequency ranges, labelled as mode 1 , 2 and 3 . Table 2 provides the frequency, magnitudes and growth rate of these modes. All these modes have comparatively large magnitudes with | | > 0.1. At the same time, these modes are the relatively darker ones in Figure 19(b) with decay rate | | < 0.03, which suggests that their effects play a role during the entire process.
For the branch with lower frequencies, mode 1 has been selected to illustrate the flow dynamics. Figure 20 shows the real part of the selected mode 1 with the isosurfaces of the pressure fluctuations at phase angle = 0 and 7 /8. At both instants, the key features of this mode from the pressure signals are the significant structures along the shock and compression waves caused by the reattachment. Additionally, the fluctuations around the shock and reattachment are three dimensional, and indicate a slight wrinkling of the shock. Comparing the modal fluctuations at these two phases, the sign switch between them describes the oscillation of the reattachment shock. Figure 21(a) provides the pressure fluctuations at the slice = 0, in which the effect that mode 1 has on shock and compression waves is more clear. Note that perturbations in the upstream turbulent boundary layer are too weak to be visible at the given levels (| / ∞ | = 0.01) of isosurfaces and in the contours.
In Figure 22, the fluctuations of the streamwise velocity component from DMD mode 1 are given. As we can see, large fluctuations are aligned with the streamwise direction with negative and positive values alternating in the spanwise direction. These longitudinal structures appear   to start within the fore part of the free shear layer and extend beyond the reattachment location. Additionally, they are mainly located in the near-wall part of the boundary layer, as shown by the streamwise velocity fluctuations in Figure 21(b). We also superimpose the modal fluctuations onto the mean flow and plot the contours of streamwise velocity in the − slice at / 0 = −2.875 in Figure 23. The high-and low-speed streaks are obvious in the contours and show very similar features as the contours of skin friction in Figure 10. Other low-frequency modes between = 0.008 ∼ 0.05 have been also examined, and they all share common structures with mode 1 . The pressure and velocity fluctuations of DMD mode 1 suggest that the low-frequency flapping motion of the shock is coupled with streamwise-elongated structures in the interaction region.
These spanwise-aligned structures are the signature of counter-rotating Görtler-like vortices, as shown by the contours of modal streamwise vorticity in Figure 24. The location and strength of  these vortex pairs are changing with the phase angles. At the given instants ( = 0 and = 5 /8), the spanwise wavelength of the vortex pair is ranging from 1.9 0 to 1.6 0 . To better visualize the unsteady motions represented by this mode, an animation of modal fluctuations with time is provided in the supplemental material (see supplementary movie 1). These time sequential snapshots of the reconstructed flow field from 1 shows the flapping motion of the reattachment shock, i.e., the displacement of the shock around its mean shock location. The counter-rotating Görtler vortices are also unsteady in terms of their strength. Additionally, figure 24 indicates that these vortices can move in both the spanwise and wall-normal directions. For mode 2 , the pressure isosurfaces in Figure 25 show high levels of fluctuations along  the reattachment shock, but the three-dimensional features are stronger compared to mode 1 . Positive and negative fluctuations are alternating in both spanwise and wall-normal directions, which represents a propagation of waves from the shear layer and outwards along the shock. The radiation of the Mach-like waves is in agreement with the results from a global linear stability analysis of an impinging shock case in a laminar regime (Guiho et al. 2016). The emission of these waves induces large disturbances along the streamwise direction in the supersonic part of the flow field. In the contours of modal spanwise-averaged pressure fluctuations in Figure 27, the radiation of the waves along the streamwise direction and shock is easier to observe. Considering the fluctuations of the streamwise velocity, shown in Figure 26, smaller longitudinal vortical structures are observed, compared to mode 1 . These vortices alternate along both the spanwise and streamwise directions, and are mainly concentrated within the boundary layer. Clearly, this mode represents the convection of the shear layer vortices and the induced Mach-like waves in the supersonic part of the flow field, which can also be seen in the contours of modal spanwise-averaged streamwise velocity in Figure 27(b). Similar observations were also reported in the two-dimensional DMD analysis of an incident shock case (Pasquariello et al. 2017). From the animation of the reconstructed flow based on this mode (see supplementary movie 2), the propagation of the Mach-like waves starting from the shock foot and the shedding of the small streamwise vortices are obvious.
The higher frequency modes, branch C, are related to the small-scale turbulent dynamics. For example for mode 3 , pressure fluctuations in Figure 28 show small highly three-dimensional arc-shaped vortices. These spanwise-aligned vortices are generated from the downstream part of the shear layer. The streamwise displacement of the fluctuations contours at different phase angles indicates the convection of the coherent vortices.
The convection behaviour of this mode is also evident from isosurfaces of the streamwise    Figure 29. The velocity fluctuations originate from the strong shear layer behind the step. It is also noticed that these fluctuations are distributed along the free shear layer and downstream boundary layer. Additionally, this mode shows less anisotropic features, compared with the other two modes. The frequency of this mode is close to the typical frequency of the turbulence considering the thicker boundary layer downstream the step. Thus, we consider this mode to be related to the convection of typical turbulent structures (see supplementary movie 3) that result from an amplification of the incoming turbulence by the separation bubble, cf. the stability analysis of Guiho et al. (2016) for an incident shock SWBLI case. We also performed a two-dimensional DMD analysis for the laminar case and similarly divided the modes into three branches (Hu et al. 2019). The branch with higher frequencies centred at 0 / ∞ ≈ 0.1 is associated with the shedding of large coherent shear vortices, which is also observed in the present turbulent case. The other two branches with lower frequencies are related to the unsteady motions of the separation bubble and the shock. In contrast to the turbulent case, we found no evidence of Görtler-like vortices in the laminar case.

Physical mechanism of low-frequency unsteadiness
The current BFS case shows unsteady motions at similar low frequencies as those usually observed for SWBLI on flat plates and on compression ramps. Compared with these cases, however, the flow topology of the present case shows significantly different features. In canonical impinging shock and ramp cases, the separation bubble is enclosed by a separation shock and reattachment shock. In contrast, the recirculation region in a BFS case is surrounded by the step expansion fan, and reattachment shock. In terms of the mean skin friction, the recirculating flow is usually less uniform downstream of the mean separation position and recovers slower in the ramp and incident shock cases (Pasquariello et al. 2017;Priebe & Martín 2012). The fluctuations of inside the separation bubble in these cases are usually related to the low-frequency unsteadiness. In the current case, however, the skin friction ( Figure 6) is relatively uniform in the upstream part of the bubble, which is caused by the 'dead-air' region close to the stationary step wall. The wall pressure is usually increasing throughout the separation bubble in the ramp and incident shock cases (Priebe & Martín 2012;Sansica et al. 2016). In the current case, the pressure drastically drops at the step and keeps a relatively steady low level in the separation bubble, which is typical for strong interactions. These differences may suggest different low-frequency features between these cases. To compare with other canonical SWBLI cases, the dimensional frequencies are rescaled based on the separation length as = / ∞ in the following discussion. The instantaneous visualizations displayed in §3.3 illustrate both relatively localized and global unsteady motions in the flow field, involving high and low speed streaks, breathing bubble and oscillating reattachment shock, as well as vortex shedding in the separated shear layer. A linear stability analysis of the mean flow shows that the most unstable global mode is mainly distributed along the dividing line ( = 0), especially around the reattachment location (Guiho et al. 2016). The RMS wall-normal velocity in Figure 8 is consistent with this observation. The spectral analysis in section 3.4 reveals that there are two kinds of low-frequency unsteadiness in the interaction region and the lower frequency around = 0.18 is associated with the coupling of separation bubble and reattachment shock wave. Furthermore, the three-dimensional DMD analysis separates the different dynamics contributing to low-frequency interaction. Apart from the unsteady separation bubble and reattachment shock, the low-frequency mode 1 also reveals the unsteady Görtler-like vortices, see Figure 24. Although these Görtler vortices are rather weak compared to other energetic dynamics such that they are hard to identify in the vortical visualization in Figure 7, the skin friction contours in Figure 10 capture the footprint of the associated high-and low-speed streaks.
These observations share qualitative similarities to the low-frequency DMD modes calculated by Priebe et al. (2016). In their ramp case, the fluctuations of the low-frequency mode clearly show shocks (mainly separation shock) and longitudinal Görtler vortices near the reattachment. In the impinging shock case of Pasquariello et al. (2017), both the visualization of streamwise vorticity and the DMD analysis of the skin friction support the finding of the Görtler vortices downstream the reattachment location. Both report that the frequency of this unsteadiness is between 0.01 < < 0.2, while the current results for a BFS correspond to a Strouhal number range of 0.03 < < 0.6, which is three times larger than the values of other canonical cases. We believe that the higher frequency in the current case is caused by the fixed separation location and confinement by the step wall. In the ramp and impinging shock cases, the recirculation region can move from both separation and reattachment sides. By comparison, the current case can only move in the reattachment part due to the limitation of the step, and it is reasonable to assume that this leads to a smaller oscillation amplitude and correspondingly to a higher frequency. This explanation is supported by the temporal evolution of the reattachment point. In Figure 15(b), the calculated minimum, mean and maximum reattachment location are / 0 = 8.3, 8.9 and 10.2, which leads to an oscillation range of about 15% . For ramp cases, oscillations of up to 70% have been reported (Priebe & Martín 2012). Moreover, the separation length is around three times the maximum separation height ( = 3ℎ, the maximum separation height equals to the step height) in our BFS case, whereas the recirculating flow regions are typically much thinner in ramp and impinging shock cases. Estruch-Samper & Chandola (2018) proposed an entrainment-recharge mechanism to associate the low-frequency unsteadiness with the shedding effects. In this theory, the Strouhal number of the low-frequency breathing can be related to the entrainment frequency by where is the length-to-thickness ratio of the shedding coherent structures; represents the percentage of the entrainment mass and is the spreading rate of the mixing layer. Huang & Estruch-Samper (2018) showed that the variations of these three parameters between different cases are small if the incoming flow conditions are close, which results in an approximate constant ≈ 0.08. The ratio of the bubble length to bubble height /ℎ and the dimensionless entrainment length ent depend on the specific geometry. In a similar entrainment and injection model by Piponniau et al. (2009), they consider the entrainment usually only occurs in the rear half of the separation bubble, i.e., the downward part of the shear layer, which leads to a dimensionless entrainment length ent ≈ 0.5 in the impinging shock and ramp cases. The geometry dependent transformation factor is ent = ent 2 /ℎ ≈ 1.5 for these canonical cases. In the current BFS case, the entrainment length ent is one, which gives the transform factor ent ≈ 3. The entrainment frequencies ent of the shear layer shedding behaviour are similar for all these cases; thus the BFS case will yield a about two times larger low than the impinging shock and ramp cases. However, this model only provides an estimate of the low frequency, and we do not expect to obtain an accurate value.
Several works in the literature have found evidence of Görtler-like vortices in SWBLI. Görtler vortices typically have a spanwise length-scale in the order of the incoming boundary layer thickness (Settles et al. 1979). The Görtler number, defined as, gives an indication on whether such vortices can form (Smits & Dussauge 2006), where is the radius of curvature of the streamline, * is the boundary layer displacement thickness and is the boundary layer momentum thickness. Figure 30 shows the curvature / and Görtler number along two streamlines of the mean flow inside the shear layer (shown in Figure 21 and Figure 27). Streamline 1 is closer to the wall and passes through the coordinate / 0 = 0 and / 0 = 0.125. Significantly large streamline curvature occurs at the separation point and around the reattachment. Correspondingly, two distinct peaks are observed at these locations. The large curvature and Görtler number at the separation is mainly caused by the sudden change of the geometry at the step edge. In a laminar flow, the critical Görtler number is around = 0.6 (marked as gray dot-dashed line) for a wide range of , above which the flow exhibits significant centrifugal instability and local Görtler vortices will emerge inside the boundary layer (Smits & Dussauge 2006). We see that the Görtler number is larger than the critical value between 7.7 / 0 18.5 and reaches its extremum = 1.2 at / 0 = 10.9, close to the reattachment. Streamline 2 is in the middle of the boundary layer and passes through the point / 0 = 0 and / 0 = 0.5625. The spatial evolution of the curvature and Görtler number has a similar trend as for the streamline 1 and the Görtler number is above the critical value for 7.0 / 0 19.3. In Figure 10, we show that the high and low speed streaks are observed within the region 6.8 / 0 17.8, which is consistent with the zones with large Görtler number. The intrinsic spanwise wavelength of these streamwise vortex pairs is reported as ≈ 2 in SWBLI systems (Schülein & Trofimov 2011;Priebe et al. 2016), which is also in agreement with the current observations = 2.0 0 ∼ 3.1 0 . The streamwise velocity field reconstructed from DMD mode 1 also displays these high and low speed streaks, see Figure 23. Although there is no general critical reported in the literature for turbulent separated flow, the high levels of provide an indication of sufficiently strong Görtler instability at the reattachment location, which provides an explanation for the low and high speed streaks observed in Figure 10 and the detected streamwiseoriented structures in DMD mode 1 .
Görtler vortices resulting from strong curvature could be unsteady in the turbulent flow, as concluded by Floryan (1991) from various experiments in low-speed turbulent flows. One of the situations proposed is that the generated streamwise-oriented vortices oscillate in the spanwise direction if they are affected by three-dimensional disturbances in the incoming flow. The unsteady streamwise vortices observed in the incident shock and ramp cases (Sun et al. 2012;Pasquariello et al. 2017;Priebe et al. 2016) both fall into this category. The mode structure shown in Figure 22 oscillates in spanwise direction. It is noticed that this proposition involves incoming disturbances, which may suggest a certain dependence on upstream flow conditions. From Figure 21(b), we can observe that weak upstream disturbances and fluctuations are part of the same DMD modes as the downstream Görtler vortices, which manifests that the observed Görtler vortices in the present study indeed have a significant correlation with upstream disturbances. For comparison, we also show the Görtler number for our laminar case (Hu et al. 2019) in Figure 31. As we can see, the curvature around the reattachment location in the laminar case is smaller than the one in the turbulent flow. As a result, has a smaller value than for the turbulent case (cf. Figure 30). Specifically, is below the critical value in the whole separation bubble, which probably is the reason that there are no Görtler vortices around the reattachment point in the laminar case (cf. Figure 12). In addition, this difference also shows how the existence of Görtler vortices is affected by upstream fluctuations: more turbulent incoming flow leads to a smaller separation length and thus stronger Görtler instability. On the other hand, the flow field around the reattachment location is more likely to reorganize and form into the spanwise-aligned vortices in the turbulent regime due to the incoming three-dimensional fluctuations.
Based on the above discussion, the following physical mechanism is proposed for the production of the low-frequency unsteadiness. The incoming turbulent flow experiences strong shear and curvature upon separation, which leads to large coherent vortical structures along the shear layer. Near the reattachment, there is significant centrifugal instability within the shear layer due to the concave streamlines. The Görtler instability, excited by incoming 3D turbulence, leads to large streamwise oriented vortices, which produce high-and low-speed streaks around the reattachment region (cf. Figure 10 and 23). These Görtler vortices are unsteady, which leads to spanwise shock wrinkling at very low frequencies, as we see from the streamwise velocity fluctuations from DMD mode 1 . Therefore, we believe that the centrifugal force and induced Görtler vortices are the main driving force of the global low-frequency unsteadiness in the turbulent case, which suggests that the low-frequency oscillation of SWBLI is inherently a three-dimensional mechanism. In the meantime, there is also notable dependence on the upstream fluctuations within the incoming turbulent boundary layer.

Conclusions
The unsteady dynamics of SWBLI over a BFS, with particular attention to the low-frequency unsteadiness, has been investigated at = 1.7 using a well-resolved LES. The mean flow field illustrates the main flow topology of SWBLI in the BFS, consisting of a centred Prandtl-Meyer expansion fan originating from the fixed separation point, a separation bubble behind the step and a reattachment shock generating from the compression waves. Different from the canonical impinging shock and compression ramp cases, the separation point is stationary and only one shock occurs in the BFS case. The instantaneous flow field shows that the unsteady behaviour is however similar to other SWBLI configurations, including the vortex shedding in the separated shear layer, as well as the breathing of the separation bubble and a flapping shock motion. The spectral analysis shows that there is a broad band of low-frequency oscillations, which we classify into two branches with the dominant frequencies centred near = 0.02 and 0.2 in the current case. The lower frequency dynamics is related to the unsteady separation bubble size, as well as the shock angle and position, while the second one connects to the shedding of shear layer vortices, which also affects the reattachment location.
Three-dimensional DMD analysis was used to reveal the characteristic mode structures that contribute to the observed unsteady behaviour. The low-frequency mode 1 provides evidence for the statistical link between the shock motion (by pressure fluctuations) and the unsteady Görtler vortices (by the streamwise velocity fluctuations) around the reattachment position. The highand low-speed streaks in the contours of in Figure 10 and the reconstructed velocity field in Figure 23 are the signature of these spanwise-aligned vortices. The medium-frequency mode 2 represents shear-layer vortices and Mach-like waves. We thus believe that the unsteady Görtler vortices around the reattachment provide the unsteady forcing that sustains the low-frequency motions of shock and separation bubble. In particular through the comparison with a laminar inflow case (Hu et al. 2019), we show that the upstream fluctuations have a notable effect on the formation and existence of the unsteady spanwise-aligned Görtler vortices. D