Bi-global stability of supersonic backward-facing step flow

Abstract

Supersonic backward-facing step (BFS) flow is numerically studied using direct numerical simulation (DNS) and global stability analysis (GSA) with a free stream Mach number of 2.16 and a Reynolds number of 7.938 × 10⁵ based on the flat-plate length L and free stream conditions. Two-dimensional BFS flow becomes unstable to three-dimensional perturbations as the step height h exceeds a certain value, while no two-dimensionally unstable mode is found. Global instability occurs with the fragmentation of the primary separation vortex downstream of the step. Two stationary modes and one oscillatory unstable mode are obtained at a supercritical ratio of L/h = 32.14, among which the two stationary modes originate from the coalescence of a pair of conjugate modes. The most unstable mode manifests itself as streamwise streaks in the reattached boundary layer, which is similar to that in shock-induced separated flow, although the flow separation mechanisms are different. Without introducing any external disturbances, the DNS captures the preferred perturbations and produces a growth rate in agreement with the GSA prediction in the linear growth stage. In the quasi-steady stage, the secondary separation vortex breaks up into several small bubbles, and the number of streamwise streaks is doubled. A low-frequency unsteadiness that may be associated with the oscillatory mode is also present.

1. Introduction

Nominally two-dimensional (2-D) surface irregularities in the form of backward-facing steps (BFSs), forward-facing steps and gaps have been a long-standing topic of interest in both theoretical research and engineering fields. Supersonic BFS flow represents a canonical case for the study of flow separation, which also occurs in many engineering applications, such as flame stabilization and supersonic mixing. Although the geometry is simple, the flow structure behind the step is not, particularly in the supersonic flow regime. Such flow is usually characterized by expansion waves, lip and reattachment shock waves, and a recirculation zone beneath a shear layer, which make supersonic BFS attractive for fundamental research of various flow physics. However, many studies have focused on incompressible BFS flow, leaving the complex flow physics in the supersonic flow regime unclear. The main motivation of the present study is to establish a primary understanding of global instability in supersonic BFS flow.

Numerous investigations have been carried out to understand subsonic BFS flow features through computational and experimental methods, which provide a relatively abundant understanding of various flow phenomena, including convective and intrinsic instabilities. The Kelvin–Helmholtz (K–H) instability, which is introduced by the shear layer generated by the separation region, is one of the most well-known convective instabilities in subsonic BFS flow, and it has mostly been studied at rather large step heights. To date, there have been considerable experimental and numerical studies that show evidence of K–H instabilities in subsonic BFS flow (Hasan 1992; Rani & Sheu 2006; Rani, Sheu & Tsai 2007; Schäfer, Breuer & Durst 2009; Schröder et al. 2013; Bolgar, Scharnowski & Kähler 2018; Eppink 2022). Kaiktsis et al. (1991) found that the BFS flow with an expansion ratio of $\varGamma $= 0.5 (i.e. the ratio of the step height to the channel height downstream) was convectively unstable to the K–H instability at Re = 525 (Re is based on half of the downstream channel height and the maximum incoming velocity on the centreline) in a large portion of the flow domain using local stability analysis based on the parallel-flow assumption. Note that their geometry is different from that in figure 4. Nowruzi, Salman Nourazar & Ghassemi (2018) adopted an energy gradient method to study the convective instability of subsonic BFS flow under different Reynolds numbers and expansion ratios, and concluded that the onset of instability originated in the separated shear layer. A prevailing viewpoint suggests that the dynamic processes observed in such subsonic BFS flow, particularly the oscillations of separated shear layers, arise from the presence of vortical structures known as K–H vortices. At higher Reynolds numbers, these K–H vortices become increasingly prominent and influential in shaping the flow patterns and dynamics within the shear layer. As these vortices are convected downstream, they induce significant alterations to the overall flow structure. Another focus of research is the step's potential to destabilize Tollmien–Schlichting (T–S) waves. As the entering T–S waves develop, especially in flows with higher Reynolds numbers, a combined instability mechanism with K–H instability downstream of the step can occur (Wang & Gaster 2005; Eppink et al. 2018; Crouch & Kosorygin 2020). Recently, Hildebrand et al. (2022) performed a linear stability analysis of a BFS flow without the confinement of the upper wall and found that the growth rate of the T–S wave significantly increased in the separation region. K–H instabilities can also be distinguished in direct numerical simulation (DNS) with a turbulent inflow in internal BFS flows (Le, Moin & Kim 1997a; Barri et al. 2010; Pont-Vílchez et al. 2019). Moreover, temporal transient growth analysis was performed by Blackburn, Barkley & Sherwin (2008) at a globally stable Reynolds number. They found that the optimal initial disturbances were concentrated around the step edge and the corresponding outcomes were periodic wavy structures.

The intrinsic instability in subsonic BFS flow has also been extensively studied. Armaly et al. (1983) reported that the two-dimensionality would be lost for a Reynolds number of approximately 400 and that the main flow instability occurs in the primary separation zone. The intrinsic instability of a 2-D flow with $\varGamma $ = 0.5 was examined by Kaiktsis, Karniadakis & Orszag (1996) using global stability analysis (GSA). The flow was shown to be globally stable to 2-D perturbations up to Re = 1875, the highest Reynolds number under consideration. Barkley, Gomes & Henderson (2002) found that the 2-D flow is unstable to three-dimensional (3-D) perturbations beyond a critical Reynolds number of 748 with $\varGamma $ = 0.5. More recently, Lanzerstorfer & Kuhlmann (2012b) completed a systematic study on the global stability of BFS flow with varying expansion ratios (from $\varGamma $ = 0.25 to $\varGamma $ = 0.975) and concluded that the instability cannot be attributed to either a centrifugal (Barkley et al. 2002) or a Taylor‒Görtler instability mechanism (Kaiktsis et al. 1996). The global stability problem with a symmetric expansion was also studied by Lanzerstorfer & Kuhlmann (2012a).

In contrast, flow instability and the transition to turbulence induced by a BFS are much less understood in the supersonic flow regime. Most studies have been conducted in fully laminar and turbulent flow regimes (Rom & Seginer 1964; Gai et al. 1989; Le, Moin & Kim 1997b; Soni, Arya & De 2017; Bolgar et al. 2018; Hu, Hickel & Van Oudheusden 2021). However, little attention has been given to the transition of supersonic BFS flow. Ginoux (1960) conducted a series of wind-tunnel experiments of supersonic BFS flow at Mach 2.16 with different flat-plate lengths and step heights while keeping the incoming boundary layer laminar. Streak-like flow structures downstream of the step were observed using a sublimation technique, which cannot be explained by upstream disturbances in the boundary layer or model imperfections, as reported in the study. Chen et al. (2013) and Zhu et al. (2015) used a nanoparticle-based planar laser scattering method to visualize the flow structure in supersonic BFS flow. The incoming boundary layer was laminar. Streak-like flow structures were also observed; however, an explanation for this phenomenon was lacking.

Only a few studies have focused on the instability of supersonic BFSs, especially the intrinsic instability. Hu, Hickel & Van Oudheusden (2019) and Hu et al. (2021) numerically studied the interaction of the incoming first mode and the K–H instability under four different inflow conditions. In the case of a laminar boundary layer with zero amplitude of the first mode, K–H vortices were observed in the shear layer formed by the separated zone downstream of the step. No streamwise streaks were found in the reattached boundary layer. The dominant frequency for the near flow reattachment is close to that captured by Schäfer et al. (2009) in the incompressible flow regime. Significantly, the incoming flow was not perturbed, indicating that an intrinsic mechanism may be responsible for the generation of vortex rolls. In fact, the separation zone caused by the sudden geometry expansion has the potential to self-excite an intrinsic instability. This has been well researched in 2-D separated flows using GSA (Theofilis, Hein & Dallmann 2000; Rodriguez & Theofilis 2010; Theofilis 2011; Taira et al. 2017). GSA is a linear stability analysis technique without making any assumptions on the spatial variation of the base flow. This feature makes GSA suitable for studying the instability of separated flows. GSA has been widely used to understand shock-induced flow separation, including supersonic/hypersonic flow over compression ramps (Hao et al. 2021; Cao et al. 2021a,b, 2022), double wedges (Sidharth et al. 2017, 2018), double cones (Hao et al. 2022), hollow cylinders/flares (Cerulus, Quintanilha & Theofilis 2021; Hao & Li 2023) and shock impingement on flat plates (Robinet 2007; Hildebrand et al. 2018; Song & Hao 2023). The abovementioned studies show that self-excited instability may occur in the separated flow if the governing parameters exceed a specific threshold (e.g. the deflection angle of the compression ramp and the strength of the incident shock wave). In addition, DNS is also used to enhance the understanding of the formation and development of the three-dimensionality, as explored by Cao et al. (2021a) and Hao et al. (2022). Both the GSA and DNS studies have revealed the existence of global instability in such separated flows.

Although various studies on BFS flow have been conducted, many details remain to be explored, particularly regarding the instability issues in supersonic flow. To the best of the authors’ knowledge, it remains unclear whether global instability exists in the separated supersonic BFS flow caused by a sudden expansion. The current research aims to investigate the intrinsic instability in supersonic BFS flow and to gain a deeper understanding of the development of the 3-D perturbations that lead to flow bifurcation using GSA and DNS. The paper is organized as follows: § 2 provides details about the numerical methodology, including DNS and GSA, the free stream flow conditions, and the geometrical parameters of the BFS. Section 3 provides the results and discussions of the base-flow, GSA and DNS. Concluding remarks are provided in § 4.

2. Numerical details

2.1. Numerical approach

The governing equations are the Navier–Stokes (N–S) equations in the following conservation form. The equations expressed in the flux-vector form are given by the equation:

(2.1)\begin{equation}\frac{{\partial \boldsymbol{U}}}{{\partial t}} + \frac{{\partial \boldsymbol{F}}}{{\partial x}} + \frac{{\partial \boldsymbol{G}}}{{\partial y}} + \frac{{\partial \boldsymbol{H}}}{{\partial z}} = \frac{{\partial {\boldsymbol{F}_v}}}{{\partial x}} + \frac{{\partial {\boldsymbol{G}_v}}}{{\partial y}} + \frac{{\partial {\boldsymbol{H}_v}}}{{\partial z}},\end{equation}

where $\boldsymbol{U = }{(\rho ,\;\rho u,\;\rho v,\;\rho w,\;\rho e)^\textrm{T}}$ is the vector of the conserved variables, ρ is the density, u, v and w are the flow velocity components in the x, y and z directions, respectively, e is the total energy per unit mass, and superscript T denotes the transpose of a matrix. Vectors F and F_v represent the inviscid and viscous fluxes in the x direction, respectively:

(2.2a,b)\begin{equation}\boldsymbol{F} = \left( {\begin{array}{*{20}{c}} {\rho u}\\ {\rho {u^2} + p}\\ {\rho uv}\\ {\rho uw}\\ {(\rho e + p)u} \end{array}} \right)\!,\quad {\boldsymbol{F}_v} = \left( {\begin{array}{*{20}{c}} 0\\ {{\tau_{xx}}}\\ {{\tau_{xy}}}\\ {{\tau_{xz}}}\\ {u{\tau_{xx}} + v{\tau_{xy}} + w{\tau_{xz}} - {q_x}} \end{array}} \right)\!,\end{equation}

where p is the pressure, and ${\tau _{ij}}$ and ${q_i}$ are the viscous shear stress tensor and heat flux vector, respectively. Vectors G, H, G_v and H_v can be expressed analogously. The equation of state for a calorically perfect gas is adopted to close the above equations. The viscosity is calculated using Sutherland's law. The specific heat ratio and Prandtl number are set to 1.4 and 0.72, respectively.

The 2-D base flow for GSA is obtained by solving the 2-D counterpart of (2.1) with an in-house multiblock parallel finite-volume solver named PHAROS (Hao, Wen & Wang 2019; Hao & Wen 2020; Hao et al. 2021). The inviscid fluxes are calculated using the modified Steger–Warming scheme (MacCormack 2014), and a monotone upstream-centred scheme for conservation laws (van Leer 1979) is implemented for reconstruction. The viscous fluxes are computed using a second-order central difference scheme. Time integration is performed by an implicit point relaxation method (Wright, Candler & Bose 1998). Far-field boundary conditions are imposed on the left and upper computational boundaries. Note that free stream conditions with no perturbations are applied to the far-field boundary. An extrapolation method is applied to the outlet boundary. The wall is assumed to be isothermal and no-slip. The 2-D simulations are initialized with the free stream conditions. The numerical solution can be considered to converge when the Euclidean norm of the density residual is reduced by approximately ten orders of magnitude. The validation of the PHAROS code for supersonic BFS flow can be found in Appendix B.

The 3-D simulation is also performed using PHAROS. The initialization method of extruding the 2-D base flow in the spanwise direction is adopted. Periodic boundary conditions are applied to the spanwise boundaries. A second-order implicit scheme is used for time discretization (Peterson 2011). The time step is set as 20 ns to ensure sufficient temporal accuracy in 3-D simulations. The other numerical schemes and boundary conditions for the DNS unsteady simulation are the same as those for the base-flow simulations.

2.2. Global stability analysis

Global stability analysis is conducted to reveal the instability of the 2-D base flow. The vector of primitive variables q = [u, v, w, T, p]^T (T and p represent the temperature and pressure, respectively) is decomposed into two parts: a base state $\bar{\boldsymbol{q}}$ and a small perturbation q′,

(2.3)\begin{equation}\boldsymbol{q}(x,y,z,t) = \bar{\boldsymbol{q}}(x,y) + \boldsymbol{q^{\prime}}(x,y,z,t),\end{equation}

where $\bar{\boldsymbol{q}}$ is assumed to be homogenous in the z direction and independent of time. The linearized N–S equations can be obtained by substituting (2.3) into the original governing equations (2.1). The perturbations are assumed to consist of infinitesimal sinusoidal disturbances:

(2.4)\begin{equation}\boldsymbol{q^{\prime}}(x,y,z,t) = \hat{\boldsymbol{q}}(x,y){\rm exp} [\textrm{i}(\beta z - \omega t)],\end{equation}

where $\hat{\boldsymbol{q}}$ is the eigenfunction, β represents the user-specified spanwise wavenumber, and ω = ω_r + iω_i contains both the growth rate ω_i (ω_i > 0, growth; ω_i < 0, decay) and the angular frequency ω_r. Here, ω_r = 0 represents a stationary mode, whereas ω_r ≠ 0 indicates an oscillatory mode. The corresponding frequency and spanwise wavelength are defined by

(2.5a,b)\begin{equation}f = \frac{{{\omega _r}}}{{2{\rm \pi} }},\quad \lambda = \frac{{2{\rm \pi} }}{\beta }.\end{equation}

Upon the substitution of q′ into the linearized N–S equations, a generalized eigenvalue problem is formed to uncover the flow instability:

(2.6)\begin{equation}\boldsymbol{\mathsf{A}}(\bar{\boldsymbol{q}};\beta )\hat{\boldsymbol{q}} = \omega \boldsymbol{\mathsf{B}}\boldsymbol{\hat{q}},\end{equation}

where $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$ are matrices associated with the spatial discretisation of the linearized N–S equations (Theofilis 2011). The generalized eigenvalue problem is solved using the implicit restarted Arnoldi method (Lehoucq, Sorensen & Yang 1998; Stewart 2002). The discretized equation (2.6) in this study is generated using an in-house GSA solver with a centred fourth-order finite-difference method. A third-order finite-difference scheme is adopted near the boundaries. Specifically, the solver adds a semi-artificial viscosity to stabilize the numerical scheme when shock waves are involved in the base flow (Hildebrand 2019). The strength of the artificial viscosity needs to be picked carefully so that it does not affect the overall physical properties of the stability problem and selectively ensures numerical accuracy. The value is of the order of 10⁻² in the present study after thorough testing and sensitivity analyses. Because the computational domain of BFS flow cannot be regarded as a combined domain, a multidomain strategy is applied in the solver according to de Vicente et al. (2006), as schematically illustrated in figure 1. The multidomain strategy extends the eigenvalue problem of the class of (2.6) to include multiple domains. In this approach, the linearized N–S equations of the disturbances are written in each of the domains connected by the interior faces. All the perturbations are assumed to be zero on the inflow boundary. The following boundary conditions are imposed on the wall surface:

(2.7a,b)\begin{equation}\hat{u} = \hat{v} = \hat{w} = \hat{T} = 0,\quad \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\hat{p} = 0,\end{equation}

where vector n represents the normal direction of the wall surface. The density can be deduced from the equation of state. On the outflow boundary, a linear extrapolation is enforced. To avoid any reflection of perturbations, sponge layers are imposed near the far-field and outflow boundaries (Mani 2012). Two validation cases for the GSA code can be found in Appendix A.

Figure 1. Geometric configuration of the backward-facing step (not to scale). The origin of the coordinate system is located at the step edge. Here h and L represent the step height and the flat plate length, respectively. Three coloured domains are used for GSA.

2.3. Geometric configuration and flow conditions

The geometrical configuration and flow conditions of the numerically studied backward-facing step in the present study are taken from the experiments conducted in the TCEA S‒1 supersonic wind tunnel by Ginoux (1960). The BFS comprises a flat plate followed by a step. The length of the flat plate is 225 mm. Five different step heights are considered, i.e. h = 5, 6, 7, 8 and 15 mm, and the case of h = 15 mm is identical to model S-1 in the experiments. The corresponding L/h ratios are 45.00, 37.50, 32.14, 28.13 and 15.00. Because the length of the flat plate is fixed, the geometrical parameters of the cases are characterized by the step height h in the following sections. The free stream has a Mach number of 2.16. The total pressure and total temperature are 29330.92 Pa and 288.15 K, respectively. The Reynolds number based on the free stream properties and flat-plate length is approximately 7.938 × 10⁵. The temperature of the wall is specified as 293 K.

A grid independence study is conducted to obtain an accurate base flow for GSA. The height of the computational inlet in front of the step is denoted by H, and the origin of the coordinate system is located at the lower edge of the step. The computational domain has an extension downstream of the step, denoted by L_d, as shown in figure 1. To demonstrate mesh convergence, three meshes are constructed with details reported in table 1. The streamwise cell numbers (N_x) are given both upstream and downstream of the step. Similarly, the wall-normal cell numbers (N_y) are shown for the computational inlet and step. Moreover, the total cell numbers are labelled by the symbol N. The computational mesh has a cluster of cells near the walls, the step, the reattachment point and the shear layer downstream of the step. The thickness of the first layer normal to the wall is approximately 1 × 10⁻⁶ m to ensure that the grid Reynolds number has an order of magnitude of one.

Table 1. Computational grid details for base flow simulations and GSA.

Figure 2 shows the distribution of the skin friction coefficient and surface pressure coefficient on the surface downstream of the step, which are obtained from three different meshes for h = 7 mm. The skin friction coefficient C_f and surface pressure coefficient C_p are defined by

(2.8a,b)\begin{equation}{C_f} = \frac{{{\tau _w}}}{{0.5{\rho _\infty }u_\infty ^2}},\quad {C_p} = \frac{{{p_w}}}{{0.5{\rho _\infty }u_\infty ^2}},\end{equation}

where ${\tau _w}$ and ${p_w}$ are the wall shear stress and pressure, respectively. The results from the medium and fine meshes show good agreement with each other. The flow region behind the step usually contains a primary vortex and a series of Moffatt vortices (Moffatt 1964). Theoretically, there are many Moffatt vortices near the corner behind the step if the grid resolution is sufficient. Herein, the length of the separation region is defined as the axial distance between the step and the reattachment point, which is the location C_f = 0 downstream of the primary vortex. In figure 2, the reattachment point is located at approximately x/h = 12.1. Quantitatively, the length of the separation zone shows a considerable difference between the results. The length on the coarse mesh is found to be approximately 6.2 % smaller than that on the fine and medium mesh. The pressure coefficient shows the same tendency. The results obtained from the medium and fine meshes almost overlap. Thus, a medium mesh resolution is adopted in the present study, which can ensure an accurate result and save computational resources.

Figure 2. Distributions of the (a) skin friction coefficient and (b) surface pressure coefficient obtained from three different meshes for h = 7 mm.

A grid independence study for GSA is also conducted using a fine mesh (330 000 cells), a medium mesh (176 000 cells) and a coarse mesh (85 500 cells). Note that the grid for computing the base flow is different from that for GSA; the flow information should be reliably interpolated on the latter grid using a proper interpolation method, which is a common procedure in processing the base flow for GSA (Theofilis, Duck & Owen 2004; Kitsios et al. 2009). Figure 3 shows the eigenvalue spectra obtained on different meshes for the case h = 7 mm corresponding to the most unstable wavenumber β = 4.5 (discussed in § 3.2). Two stationary modes and one oscillatory unstable mode are well captured by the medium and fine mesh. The relative error of the growth rate of the most unstable mode between fine and medium mesh is approximately 0.84 %. As a result, the medium mesh is sufficient for GSA.

Figure 3. Comparison of eigenvalue spectra obtained on different meshes for the case h = 7 mm corresponding to the most unstable wavenumber β = 4.5. Square, fine mesh; diamond, medium; triangle, coarse mesh.

3. Results

3.1. Base flow features

In this section, the base flow fields obtained from different cases are discussed. The overall flow structure of the supersonic BFS flow can be observed in figure 4. A leading-edge shock wave is generated at the lip of the flat plate owing to the viscous interaction (Anderson 2000). Afterwards, the supersonic flow encounters the step and an expansion fan is generated at the corner, resulting in a pressure drop behind the step. The flow separates behind the step, forming a recirculation zone and a lip shock. In the recirculation zone, there is a primary vortex located upstream of the reattachment point and a secondary vortex near the origin, rotating in opposite directions. The two vortices are separated by the secondary separation point on the wall, which is the position at which the primary vortex separates from the surface downstream of the step. Additionally, the flow velocity in the recirculation zone is relatively smaller than that of the main flow. A shear layer develops between the free stream and the recirculation zone. As the shear layer approaches the downstream wall, the flow is compressed to generate the recompression shock wave. After the reattachment point, the shear layer promotes the redeveloped boundary layer. Overall, the behaviour of the flow in BFS can be described as follows: the external flow changes direction by turning downward at a right angle due to the expansion waves. This leads to the formation of a recirculation zone with several vortices inside. Then, the flow almost returns to its original direction parallel to the free stream flow due to the recompression shock wave.

Figure 4. Schematic of the flow structure over a backward-facing step.

Contours of the Mach number and pressure near the primary vortex superimposed with the streamlines inside the recirculation zone and the isoline of u = 0 are shown in figure 5 for different cases. The isoline of u = 0 denoted by the black dashed line is adopted to characterize the recirculating zone. Notably, the range of the pressure contour in figure 5(f) differs from the others to clearly show the pressure variations. The incoming flow expands at the step corner due to the sudden geometry expansion. For all cases, the expansion angles are approximately identical; thus, there are slight differences in the Mach number contours. There are two vortices inside the recirculation zone. A sequence of Moffatt vortices can be found if sufficient grid cells are placed near the concave corner of the step. In the present study, more than one Moffatt vortex can be distinguished under the adopted grid resolution. Based on Theofilis et al. (2004) and the validation case presented in Appendix A, it has been established that the presence of small Moffat vortices in the lid-driven flow does not significantly impact the flow stability. Consequently, we will not discuss the behaviour of these small vortices in the following sections, apart from the primary and secondary vortices. Exploring the detailed behaviour of these additional vortices is beyond the scope of this paper.

Figure 5. Contours of the Mach number and pressure for supersonic BFS flow with step heights ranging from h = 5 to 15 mm. The left column represents Mach number contours superimposed with the u = 0 isoline (black dashed line) in the recirculation zone with streamlines (green solid line), and the pressure contour near the primary vortex with streamlines (green solid line) is shown in the right column.

Generally, the non-dimensional length of the recirculation zone based on h decreases with increasing step height. As the secondary vortex grows, the primary vortex becomes slenderer and has a more preferential orientation along the shear layer. The primary vortex exhibits a smooth oval shape in the case of h = 5 mm, which then becomes distorted with increasing step height. At the step height of h = 6 mm, the core of the primary vortex cannot maintain a typical oval shape. One side of the core in the longitudinal direction becomes accentuated and constricts. As the step height reaches h = 7 mm, the core distorts completely, and by h = 8 mm, it fragments into two distinct vortices. Furthermore, under certain conditions, three vortices emerge within the primary vortex core. This fragmentation behaviour of the primary vortex is commonly observed in supersonic compression corner flows as ramp angles increase, as documented in several studies (Shvedchenko 2009; Egorov, Neiland & Shredchenko 2011; Gai & Khraibut 2019; Hao et al. 2021). Within the range of cases considered, specifically from the step height of h = 5 to 15 mm, it is essential to note that the flow remains two-dimensionally stable, with no indications of self-induced unsteadiness. However, if the step height continues to increase, it is anticipated that the flow in 2-D laminar simulations for supersonic BFS will become unstable. Additionally, when examining the streamlines close to the wall, especially near the x/h = 6 region, a slight bump can be detected at h = 7–8 mm. This bump indicates the potential formation of a minor vortex in that zone. As for the secondary vortex, its occupied space expands as the step height increases. Unlike the primary vortex, the secondary-vortex fragmentation only occurs at a greater step height, specifically at h = 15 mm in the present case. Additionally, this fragmentation behaviour tends to occur along the vortex side, rather than only within the vortex core.

Furthermore, examining the pressure distribution near the primary vortex reveals how small vortices are generated. At h = 6 mm, the pressure distribution is approximately uniform around the core of the primary vortex, which indicates that the intensity of the vortex is relatively low under the influence of the centrifugal force balanced by the pressure. An elongated low-pressure region perpendicular to the streamwise direction forms and crosses the core of the primary vortex at h = 6 and 7 mm. This low-pressure region also extends beyond the separation region, manifesting as a weak expansion wave. As h is further increased, an additional low-pressure area forms upstream of the core, leading to the observation of three distinct low-pressure regions when the step height h equals 15 mm. This observation suggests that the pressure differential creates an internal force that influences the evolution of the primary vortex's core, enhancing its development. This phenomenon is strengthened with increasing step height. Although not shown here, the pressure difference also exists in the core of the secondary vortex with a much lower intensity. As seen later, the occurrence of global instability is linked with the evolution of the primary separation vortex with increasing step height.

Figure 6 shows the distributions of the surface skin coefficient for different step heights. The open circles in figure 6 represent the reattachment points. The enlarged view near x/h = 6 is shown in figure 6(b). For all cases, the overall tendency of the skin friction coefficient is to decrease towards a global minimum value at x/h ≈ 8 first, followed by a continuous increase. Downstream of the reattachment point, the skin friction coefficient achieves the global maximum value and then remains almost unchanged at ${C_f} = 5 \times {10^{ - 4}}$. At x/h ≈ 7, C_f fluctuates if the step height is larger than 6 mm, which corresponds to the observation in figure 5 where the streamlines in this region become humped, indicating that a near-wall vortex will be generated. In particular, the near-wall vortex can be seen at h = 15 mm. Figure 7 shows the distribution of the length of the recirculation zone with the step height. It is confirmed that the non-dimensional length of the recirculation region decreases with increasing step height. However, the dimensional length increases linearly with the step height, as shown by the black dash-dotted line in figure 7.

Figure 6. (a) Distributions of the surface skin coefficient for different step heights. (b) Enlarged view at x/h ≈ 6. Open circle, reattachment points. The horizontal dashed dotted line indicates zero skin friction.

Figure 7. Distributions of the length of the recirculation zone with step height h. Red solid line, nondimensional length; blue solid line, dimensional length; black dash-dotted line, linear curve fitting.

Figure 8 shows the distributions of the surface pressure coefficient for different step heights. Figure 8(b) is the enlarged view near x/h = 6. The pressure in the plateau region located at x/h ≤ 6 decreases with increasing step height. At h = 5 mm, the plateau region is immediately followed by a pressure rise that reaches a constant value of C_p = 0.31 at x/h > 13. This is slightly larger than the free stream pressure under the influence of the weak recompression shock. However, for h ≥ 6 mm, the pressure first decreases to its global minimum and then increases to its constant value, which can be confirmed by figure 8(b). In fact, the global minimum of surface pressure corresponds to the footprint of the low-pressure region observed in figure 5. Thus, the pressure tendency for h ≥ 6 mm poses a local adverse pressure gradient to the reverse flow, which promotes the separation of the bottom flow and the generation of the near-wall vortex.

Figure 8. (a) Distributions of surface pressure coefficients at different step heights. (b) Enlarged view at x/h ≈ 6. Open circles, reattachment points.

3.2. Occurrence of global instability

Based on the analysis in § 3.1, the flow structure of the primary vortex at h = 6 mm becomes distorted according to the streamline around the core of the primary vortex. Thus, GSA is performed for the base flows at h = 5, 6 and 7 mm over a wide range of spanwise wavelengths to determine the global stability of supersonic BFS flow. The determination of the critical step height that leads to marginal instability can be found in Appendix C.

Figure 9 shows the growth rates of the two most unstable modes labelled modes 1 and 2 as a function of spanwise wavelength at different step heights. The right column of figure 9 is the enlarged view between λ/h = 0 and 2. The square and triangle symbols represent the stationary (ω_r = 0) and oscillatory modes (ω_r ≠ 0), respectively. At h = 5 mm, every eigenvalue has a negative growth rate, which means that the system is globally stable to 3-D perturbations with different spanwise wavelengths. At h = 6 mm, the supersonic BFS flow becomes unstable, with the largest growth rate occurring at λ/h = 2.094. Mode 1 is oscillatory and becomes dominant over a large range of spanwise wavelengths. However, mode 1 and its conjugate mode merge into two stationary modes at approximately λ/h = 0.6. In terms of the two stationary modes, mode 1 with a large growth rate is labelled mode 1L, while mode 1 with a small growth rate is labelled mode 1S. As the step height is further increased, the mode coalescence occurs at a larger spanwise wavelength of λ/h = 3.142. One branch of the real eigenvalues becomes the dominant mode. The maximum growth rate increases significantly, shifting to a smaller wavelength at λ/h = 1.396. Similar mode coalescence was reported in hypersonic compression corner flow (Cao et al. 2021b, 2022; Hao et al. 2021); however, no long-wavelength unstable mode can be found in the BFS flow. The other branch has a smaller growth rate; it is still unstable at λ/h = 1.396, but has a marginal impact on the instability of the system. The GSA result reveals that the critical step height for supersonic BFS at the Reynolds number considered in the present study is approximately 6 mm. The case of h = 6 mm is weakly unstable and dominated by an oscillatory mode, while the case of h = 7 mm is strongly unstable with a significantly increased growth rate and dominated by a stationary mode.

Figure 9. Growth rates of modes 1 and 2 as a function of spanwise wavelengths for different step heights: (a) h = 5 mm; (b) h = 6 mm and (c) h = 7 mm. The enlarged views around λ/h ≈ 1 are shown in the right column. Square symbols, stationary modes; triangle symbols, oscillatory modes. The horizontal line indicates the zero-growth rate; the vertical line indicates mode 1's maximum value position.

Figure 10 shows the real parts of the spanwise velocity perturbations $\hat{w}/{u_\infty }$ at the least stable wavelength for different step heights, and the corresponding eigenvalue spectra are shown in the right column of figure 10. For all cases, the spanwise velocity perturbations are mostly confined to the recirculation zone. At h = 5 mm, the least stable mode presents as a pattern of alternating positive and negative components located in the region of x/h ≤ 6. With increasing step height, the amplitude of the spanwise perturbation located at x/h ≤ 6 weakens and eventually disappears. A new pair of positive and negative components of the spanwise velocity perturbations occurs at approximately 6 ≤ x/h ≤ 8, where the skin coefficient has a local minimum value, as shown in figure 6. The step height of h = 6 mm is the critical case, which has an oscillatory unstable mode. Comparing figures 10(a) and 10(b) shows that the oscillation is probably associated with the spanwise velocity perturbations in the primary vortex at x/h ≤ 6 and that the instability is featured by the spanwise velocity perturbations at 6 ≤ x/h ≤ 8. This speculation can be proven further by the case of h = 7 mm, where the perturbations contributing to oscillation diminish at x/h ≤ 6. Meanwhile, the perturbations located at 6 ≤ x/h ≤ 8 are highly distorted as more than one pair of positive and negative velocity perturbation components. It can also be observed that the perturbation is transmitted to the free stream as oblique waves.

Figure 10. Contours of the real part of spanwise velocity perturbations (left column) for the least stable spanwise wavelength and the corresponding eigenvalue spectra (right column): (a) h = 5 mm, λ/h = 4.189; (b) h = 6 mm, λ/h = 2.094; and (c) h = 7 mm, λ/h = 1.396. Red, positive velocity perturbations; blue, negative velocity perturbations. Black dashed line, isoline of u = 0.

Figure 11 shows the real parts of the spanwise velocity perturbations of mode 1 from the lower branch (i.e. mode 1S) and mode 2 at h = 7 mm. Mode 1S exhibits the same behaviour as mode 1L, characterized by small perturbations at 6 ≤ x/h ≤ 8. Mode 2, which corresponds to an oscillatory mode, has a perturbation structure within the primary vortex similar to that of mode 1 (oscillatory) at h = 6 mm. However, the variations in the perturbations in the primary vortex of mode 2 at h = 7 mm are stronger, which appears as more than one pair of positive and negative portions. From the above discussion, it can be inferred that the global instability of supersonic BFS flow is mainly attributed to the distortion near the upstream of the reattachment point, i.e. fragmentation of the core of the primary vortex.

Figure 11. Contour of the real part of spanwise velocity perturbations at h = 7 mm for (a) mode 1S and (b) mode 2 at λ/h = 1.396. Black dashed line, isoline of u = 0.

To understand the instability mechanism of supersonic BFS flow, the GSA result is imposed onto the base flow to reveal the evolution of the most unstable mode at h = 7 mm, as shown in figure 12. The compressible stream function ψ in the recirculation zone is computed using the definition of $\rho \boldsymbol{u} = \boldsymbol{\nabla } \times \boldsymbol{\psi },\;\boldsymbol{\psi } = [0,0,\psi ]$ (Hildebrand et al. 2018) to characterize the primary vortex. The stream function can be computed through path integration of the u and v components, under the condition that the spanwise component of the velocity (w) is much smaller than the other two velocity components. Additionally, the influence of the density can be neglected since it remains nearly constant within the separation zone. By computing the stream function in each xy plane, a 3-D contour of $\psi $ can be formed. As the global mode is stationary, the isosurface of the stream function ψ = 0, denoted by the grey isosurface in figure 12, can indicate the streamwise extent of the recirculation zone.

Figure 12. Main flow structures of BFS flow reconstructed from 2-D base flow and GSA results at β = 4.5 for the case of h = 7 mm. In all figures, the grey isosurface represents the boundary of the primary vortex. (a) Blue and red isosurfaces represent the positive and negative streamwise vortices, respectively, including contours of the density gradient magnitude on the left; (b) spanwise array of the isosurfaces of the positive and negative streamwise perturbations; and (c) pressure perturbations superimposed on the grey isosurface.

Under the influence of the velocity perturbations, the straight reattachment line is distorted into a zigzag shape. The main streamwise vortices underneath the primary vortex are located in the region of 6 ≤ x/h ≤ 8, which presents as a couple of three counterrotating vortices that contribute to the distortion of the reattachment line. Because the base flow has zero spanwise velocity, the formation of the counterrotating vortices within the primary vortex are mainly induced by spanwise velocity perturbations associated with the 3-D global mode. The counterrotation of these vortices plays a role in transferring energy. These vortices attract high-momentum fluid outside the recirculation zone towards the wall and push low-momentum fluid inside the recirculation zone upwards into the flow. The couple of counterrotating vortices appear periodically as a spanwise array in the spanwise direction, further changing the base flow pattern. Furthermore, these vortices alter the smoothness of the isosurface of ψ = 0, and an uneven isosurface can be observed downstream where there is the root of the reattachment shock. The same features are shown in figure 12(c), and the velocity vector and pressure are superimposed on the isosurface of ψ = 0. The staggered pressure distributions in both the spanwise and streamwise directions can be seen clearly, which supports the spanwise momentum exchange. From another point of view, the distribution of the pressure shows evidence of the corrugation of the isosurface of ψ = 0 located at the foot of the reattachment shock. Based on the analysis of the reconstructed flow, it can be concluded that the spanwise velocity perturbations can amplify to a significant degree, leading to the development of streamwise vortices. These vortices alter the flow pattern in the spanwise direction and effectively transport the high-momentum flow into the separation zone. Furthermore, these vortices and pressure perturbations support the corrugations of the separation zone, which also contribute to the instability of the system. Nonetheless, further investigation is still needed to determine the primary source responsible for exciting the initial perturbations that ultimately induce the instability of the flow system.

To further understand the flow mechanisms that drive the induced global instability by the separation zone, a kinetic energy budget analysis is performed for the most unstable mode at h = 7 mm. In fact, the kinetic energy contributes approximately 90 % of the Chu energy (Chu 1965). Chu energy represents the total energy in the disturbance, which characterizes the intensity of the disturbances in the sense,

(3.1)\begin{equation}E = \int_\nu {\left( {\frac{1}{2}\bar{\rho }{u^{\prime}_i}{u^{\prime}_i} + \frac{1}{2}\frac{{{{\bar{a}}^2}}}{\gamma }\frac{{{{p^{\prime}}^2}}}{{\bar{\rho }}} + \frac{1}{2}\frac{{\bar{\rho }{C_v}{{T^{\prime}}^2}}}{{\bar{T}}}} \right)} \,\textrm{d}\varOmega .\end{equation}

Thus, only the kinetic disturbance energy is considered here, which is expected to provide insight into the physical instability mechanism. The linearized governing equation of velocity perturbations is given by

(3.2)\begin{equation}\frac{{\partial {u^{\prime}_i}}}{{\partial t}} + {\bar{u}_j}\frac{{\partial {u^{\prime}_i}}}{{\partial {x_j}}} =- \frac{{(\rho {u_j})^{\prime}}}{{\bar{\rho }}}\frac{{\partial {{\bar{u}}_i}}}{{\partial {x_j}}} - \frac{1}{{\bar{\rho }}}\frac{{\partial p^{\prime}}}{{\partial {x_i}}} + \frac{1}{{\bar{\rho }}}\frac{{\partial {{\tau^{\prime}_{ij}}}}}{{\partial {x_j}}},\end{equation}

where indices i = 1, 2 and 3 represent the x, y and z directions, respectively. The specific kinetic disturbance energy is defined by

(3.3)\begin{equation}{e^{\prime}_k} = {\textstyle{1 \over 2}}u^{\prime}u^{\prime {\dagger} }_i.\end{equation}

Generally, the kinetic disturbance energy ${e^{\prime}_k}$ can be obtained by multiplying (3.2) and ${u^{\prime {\dagger}}_i} $; then, one can obtain the following equation:

(3.4) \begin{equation} \frac{{\partial {e^{\prime}_k}}}{{\partial t}} = {u^{\prime}_i}^{\dagger} \frac{{\partial {u^{\prime}_i}}}{{\partial t}} = Re \left( {\underbrace{{ - {u^{\prime}_i}^{\dagger} {{\bar{u}}_j}\frac{{\partial {u^{\prime}_i}}}{{\partial {x_j}}}}}_{{\textit{Convection}}}\underbrace{{ - {u^{\prime}_i}^{\dagger} {u^{\prime}_j}\frac{{\partial {{\bar{u}}_i}}}{{\partial {x_j}}} - \frac{{\rho^{\prime}}}{{\bar{\rho }}}{u^{\prime}_i}^{\dagger} {{\bar{u}}_j}\frac{{\partial {{\bar{u}}_i}}}{{\partial {x_j}}}}}_{{\textit{Production}}}\underbrace{{ - \frac{{{u^{\prime}_i}^{\dagger} }}{{\bar{\rho }}}\frac{{\partial p^{\prime}}}{{\partial {x_i}}}}}_{{\textit{Transfer}}}\underbrace{{ + \frac{{{u^{\prime}_i}^{\dagger} }}{{\bar{\rho }}}\frac{{\partial {{\tau^{\prime}_{ij}}}}}{{\partial {x_j}}}}}_{{\textit{Viscous}}}} \right)\!,\end{equation}

where $^{\dagger} $ and Re denote the complex conjugate and the real part of a complex number, respectively. Recall that the perturbations can be written in terms of eigenmodes, as given by (2.4). The right-hand side of (3.4) represents the contributions to the time rate of change in the kinetic disturbance energy from various mechanisms. Note that the base flow is unstable if the left-hand side of (3.4) is positive and vice versa. The first term on the right-hand side of the equation represents the convection of the kinetic disturbance energy by the base flow. The second term represents the production of the kinetic disturbance energy extracted from the base flow. The third and last terms represent the work done by the perturbed pressure and viscous forces, respectively. Among the four terms, the production term contributes most to the instability of the base flow, which is consistent with the findings of Hao et al. (2021). Figure 13 shows the spatial distribution of the production term for the most unstable mode and mode 2 for h = 7 mm at the spanwise wavelength λ/h = 1.396. Clearly, the production term transfer from the base flow to the mode is mostly near the core region of the primary vortex, which confirms that the main instability arises from the separation zone. Based on the distribution of the production term, it is likely that the instability of the present supersonic BFS flow is dominated by an elliptic instability mechanism. The elliptic instability is commonly associated with viscous vortices that have elliptical streamlines (Pierrehumbert 1986, Sipp & Jacquin 1998, Kerswell (2002). Kuhlmann, Wanschura & Rath (1998) expanded this concept, proposing that the elliptical instability can also apply to strained vortices with closely aligned streamlines, referred to as a generic elliptical vortex. At the step height of h = 7 mm, the flow structure exhibits a topology similar to an elliptical vortex. Additional support for the elliptic instability mechanism comes from the observed energy transfer near the central elliptic stagnation point. In figure 13, the disturbance energy transferred from the basic flow predominantly concentrates at the centre of the vortex, particularly in the dominant mode 1L. These findings support the existence of elliptical instability.

Figure 13. Spatial distributions of the production terms in the kinetic disturbance energy equations obtained from different modes with streamlines superimposed at h = 7 mm with spanwise wavelength λ/h = 1.396: (a) mode 1L and (b) mode 2.

In summary, the spatial distribution of energy transfer, which acts as a destabilizing factor, closely resembles the pure elliptical instability flow and the subsonic BFS flow studied by Lanzerstorfer & Kuhlmann (2012b). Furthermore, the perturbations of mode 1 are primarily concentrated near the elliptical vortex, leading to the breakup of quasi-elliptical regions in the streamlines. These two distinct features, which serve as the characteristic signatures of an elliptic instability, are clearly evident in figures 10(c) and 13. Consequently, we propose that the dominant mechanism driving this instability is likely the elliptic instability, which is characterized by a peak in energy production associated with the vortex core and a highly perturbed flow aligned with the principal direction of the vortex. Moreover, the production term obtained from mode 1L is much larger than that from mode 2, which also indicates that mode 1 is the dominant mode to destabilize the base flow.

3.3. DNS of supersonic BFS flow

In this section, DNS is performed to reveal the 3-D flow structure and verify the GSA result. The case of h = 7 mm is considered here, which has two stationary and oscillatory modes. Since the most unstable mode at h = 7 mm occurs at the spanwise wavelength of $\lambda /h = \textrm{1}\textrm{.396}$ (corresponding to β = 4.5), four times the wavelengths are extended in the spanwise directions for the computational region. The distributions of grid cells in the streamwise and wall-normal directions remain unchanged as the adopted medium grid resolution, which is used for computing base flow. Seventy cells are placed uniformly for every spanwise wavelength, resulting in a total of 280 cells in the spanwise direction. The total number of grid cells for DNS is approximately 90 million. The 3-D simulation starts from tu_∞/h = 0 to 1510 (i.e. the physical time is 0–20 ms) to obtain a quasi-steady state to resolve the development of three-dimensionality and low-frequency unsteadiness. In the present study, no extra perturbations are imposed in the simulation. It is expected that the growth of instability waves can arise from the exceedingly small perturbations provided by the numerical round-off errors to depict the intrinsic instability predicted by GSA.

To characterize the development of 3-D perturbations with respect to the base flow in terms of the specific spanwise wavelength, the temporal evolution of the spanwise velocity perturbations is extracted for comparison. As illustrated by Cao et al. (2021b, 2022) and Hao et al. (2021), the averaged spanwise velocity magnitude normalized by the free stream velocity in a streamwise normal plane is calculated using the following equation:

(3.5)\begin{equation}{\sigma _w} = \sqrt {\frac{1}{{{N_y}{N_z}}}\sum\limits_{j = 1}^{{N_y}} {\sum\limits_{k = 1}^{{N_z}} {\left( {\frac{w}{{{u_\infty }}}} \right)_{j,k}^2} } } ,\end{equation}

where N_y and N_z are the numbers of grid cells in the wall-normal and spanwise directions, respectively. Figure 14 shows the temporal evolution of the averaged spanwise velocity magnitude in the z–y plane located at x/h = 7.14. This plane is located in the region where the spanwise perturbations are confined; this region represents where the main perturbations occur, as obtained from the GSA result, as shown in figure 10(c). After an initial adjustment owing to the numerical method, the exponential growth rate obtained from DNS shows good agreement with that of the most unstable mode predicted by GSA until tu_∞/h = 405, which verifies the fidelity of both GSA and DNS. The initial adjustment could be attributed to slight differences in the numerical methods employed between the 2-D and 3-D simulations. The inherent complexes, additional dimensions and applications of periodic boundary conditions in 3-D simulations can lead to slight discrepancies compared with the flow field extruded from 2-D base flow during the initial transient, which may manifest as the initial adjustment. However, it is well noted that the adjustments can be principally excluded by adapting the numerical method (Hildebrand et al. 2018), but more computational time will be needed to achieve the linear growth stage. After the exponential growth of σ_w, nonlinear saturation occurs until tu_∞/h ≈ 900, and then a quasi-steady state is established. This comparison provides strong evidence that GSA can accurately predict the most unstable mode induced by the global instability of the flow and proves the significant effectiveness of both GSA and DNS.

Figure 14. Temporal evolution of the average absolute spanwise velocity perturbations at x/h = 7.14 obtained from the DNS result compared with the growth rate of the most unstable mode predicted by GSA. Solid line, DNS; dashed line, GSA.

Figure 15(a) shows the contour of the spanwise velocity in the x–y plane of z/h = 0.7143 at tu_∞/h = 362.5 in the linear growth stage. The main flow structure is nearly identical to the spanwise velocity perturbations of mode 1L in the case of h = 7 mm predicted by GSA, as shown in figure 10(c), which also confirms that mode 1 dominates the linear growth stage of the 3-D perturbations. Figures 15(b) and 15(c) show the contours of the spanwise velocity at two y–z planes of x/h = 6.8281 and 7.6814, respectively. There is a slight attenuation of the spanwise velocity in the middle of the plane. This phenomenon is likely a result of perturbations that tend to arise from the periodic boundary conditions applied in the spanwise direction. However, this attenuation does not have a significant impact on the growth of the instability wave, as evidenced by figure 14. The perturbations have a periodic distribution with a wavelength of approximately λ/h = 1.3. For a quantitative comparison, the spatial characteristics of the perturbations are provided by the (spatial) power spectral density (PSD) of C_f along the streamline direction, as shown in figure 16. From the PSD of C_f, the peaks at the spanwise wavelength λ/h ≈ 1.3 can be extracted in the region of x ≈ 7, which is more energetic and corresponds well to the region with large spanwise velocity perturbations predicted by GSA. The wavelength λ/h obtained from the PSD is approximately 1.35, which is slightly smaller than the predetermined wavelength λ/h = 1.3916 from the GSA. A closer examination of figure 9 gives an explanation that the growth rates and the structures of mode 1 in the wavelength range of 1.23~1.39 are nearly identical. This observation suggests that the signals detected by DNS may be subject to blurring, resulting in a slight shift in wavelength.

Figure 15. Contours of the spanwise velocity at tu_∞/h = 362.5: (a) x–y plane at z/h = 0.7143; (b) y–z plane at x/h = 6.8281, denoted by A in panel (a); and (c) y–z plane at x/h = 7.6814, denoted by B in panel (a).

Figure 16. Power spectral density of the skin friction coefficient on the wall downstream of the step obtained from DNS at tu_∞/h = 362.5.

Figure 17 shows the evolution of the spanwise velocities in the later linear stage in an x–y plane of z/h = 0.7143. At tu_∞/h = 407.8 and 438.0, the spatial flow structure is still almost identical to that of mode 1 captured by the GSA. The spanwise velocity magnitude increases with time and occupies a larger region. At tu_∞/h = 438.0, the perturbations begin to propagate upstream within the primary vortex, resulting in the formation of streak-shaped structures aligned in the streamwise direction. After tu_∞/h = 468.2, the streak-shaped structures are strongly distorted and slowly break up into smaller structures. At tu_∞/h = 528.7, the flow structure looks like a combination of unstable modes captured by GSA (figures 10 and 11), which hints at the emergence of mode 2 and mode 1S. The non-dimensional frequency of mode 2 from the GSA result is approximately 0.005, which is of the same order as the low frequency (0.009) reported in the large eddy simulation (LES) conducted by Riley, Ranjan & Gaitonde (2021). Although the Reynolds number and Mach number in the study of Riley et al. (2021) are higher than those in the present study, it can be speculated that mode 2 can persist in supersonic BFS flow.

Figure 17. Temporal evolution of the spanwise velocities at the x–y plane of z/h = 0.7143 obtained from different times: (a) tu_∞/h = 377.6; (b) tu_∞/h = 407.8; (c) tu_∞/h = 438.0; (d) tu_∞/h = 468.2; (e) tu_∞/h = 498.5 and (f) tu_∞/h = 528.7.

Owing to the existence of global unstable modes, the saturated flow obtained from DNS exhibits a complicated structure in both spatial and temporal space. To further understand the influence of the perturbation on the flow structure, the C_f distributions downstream of the step at different times are shown in figure 18. Isolines of C_f = 0 are used to indicate the separation regions. Note that the vertical solid line at x/h ≈ 3.5 represents the secondary separation line (also see figure 4), which separates the primary–secondary vortices on the wall in the streamwise direction. The second solid line at x/h ≈ 12 is the reattachment line. In the initial stage, the two vertical lines present 2-D characteristics. At the later stage of the linear growth period, the reattachment and secondary separation lines become slightly distorted, but the separation zone retains the original pattern. With increasing time, the reattachment line twists further with a zigzag pattern, and similar features were reported by Cao et al. (2021b, 2022) in their DNS of hypersonic compression corner flows and by Hao et al. (2022) in their DNS of hypersonic double cone flows. Meanwhile, the secondary separation line becomes highly disordered and breaks up into several small separation zones, which are not observed in 2-D simulations. After tu_∞/h = 1057.3, the reattachment line presents a harmonic distribution characteristic, exhibiting a wavelength of λ/h ≈ 0.68, which is consistent with the second harmonic of the wavelength of mode 1. Additionally, the line oscillates at a relatively low frequency. Meanwhile, the secondary separation line is fully distorted and cannot be intuitively distinguished. More small bubbles with irregular shapes can be observed and are distributed randomly in the original secondary separation region. The peak of C_f is much higher than that in the linear growth stage, and the streamwise streaky patterns persist downstream of the reattachment line. The streak-like flow structures attest to the phenomena observed by Ginoux (1960) and Zhu et al. (2015). Thus, it can be concluded from the GSA and DNS results that intrinsic instability is responsible for the generation of streak-like flow structures. Figure 19 shows the power spectral density of the averaged spanwise velocity (σ_w) at x/h = 7.14 from th/u_∞ = 377.6 to 1510.5. Based on the spectrum, the quasi-steady BFS flow exhibits a broadband feature covering the frequency of mode 2. The low-frequency broadband characteristics account for the oscillation of C_f, as shown in figure 18.

Figure 18. Contours of the skin friction coefficient at (a) tu_∞/h = 113.3, (b) tu_∞/h = 422.9, (c) tu_∞/h =  604.2, (d) tu_∞/h = 1057.3, (e) tu_∞/h = 1283.9 and (f) tu_∞/h = 1495.4. Black solid lines, isolines of C_f = 0.

Figure 19. Power spectral density of the averaged spanwise velocity at x/h = 7.14 from tu_∞/h = 377.6 to 1510.5.

To closely examine the 3-D flow structures in the separation zone, figure 20 shows the contours of the streamwise velocity u/u_∞ (<0.99) on several y–z planes (seven equidistant distributed planes at 6 ≤ x/h ≤ 8) and one x–z plane (at y/h = 0.1) at tu_∞/h = 422.9 and tu_∞/h = 1283.9, respectively. The isolines of u/u_∞ (<0) are superimposed in the equidistantly distributed streamwise planes. The contours of the spanwise velocity are shown on the left. To facilitate processing, the isosurfaces of u = 0 are plotted to characterize the separation zones. At tu_∞/h = 422.9, the spanwise perturbations gathered at approximately 6 ≤ x/h ≤ 8 and the primary vortex with a zigzag reattachment line remained smooth. The flow structures are nearly identical to those reconstructed using the global mode and the base flow, as shown in figure 12(a). In the quasi-steady state, the spanwise velocity perturbations propagate upstream and pervade the whole separation zone. The primary vortex is severely corrugated. Several small bubbles form with irregular shapes and are randomly distributed in the secondary vortex, which reinforces the observations based on the skin friction coefficient plotted in figure 18. In addition, streamwise streaks can be observed and do not breakdown downstream. Whether the intrinsic instability is able to trigger a laminar‒turbulent transition at h = 15 mm will be the focus of a future study. That, again, would imply that the streak-like flow structures can be attributed to the intrinsic instability of supersonic BFS flow, which is two-dimensionally stable and destabilized under the influence of 3-D perturbations with a specific spanwise wavelength.

Figure 20. Contours of the streamwise velocity shown on a wall-parallel plane at y/h = 0.1 close to the wall and seven equidistantly distributed streamwise planes superimposed with isolines of streamwise velocity in the range of 0 ≤ x/h ≤ 12 at different times. Contours of the spanwise velocity at the z = 0 plane are shown in panels (c) and (d), respectively. Grey isosurfaces of u = 0 are plotted to represent the separation zones. (a) and (c) tu_∞/h = 422.9; (b) and (d) tu_∞/h = 1283.9. (The cut-off levels are u/u_∞ > 0.99 and u/u_∞ > 0 for contours and isolines of streamwise velocity, respectively.)

4. Conclusion

In the present study, DNS and GSA are performed for a supersonic backward-facing step flow where the step heights vary from 5 to 15 mm. The free stream Mach number and Reynolds number based on the flat-plate length are 2.16 and 7.94 × 10⁵, respectively. Owing to the sudden expansion downstream of the step, a large recirculation zone forms, which mainly includes a primary and a secondary zone. More Moffatt vortices can also be observed with sufficient grid resolution near the origin of the step corner. Importantly, the flows for all step heights considered in the present study are stable in a 2-D state; that is, no two-dimensionally unstable mode is found.

GSA shows that the supersonic BFS flow is unstable to 3-D perturbations when the step height h exceeds 5.4 mm in the present study. As the step height approaches 6 mm, a pair of conjugate modes are detected. Then, the pair of conjugate modes merges into two stationary unstable modes at h = 7 mm, and another oscillatory unstable mode is revealed. The spanwise velocity perturbations of the most unstable mode are mainly confined to a small region near the primary vortex core that is located close to the reattachment point, whereas the velocity perturbations in the unstable oscillatory mode almost fill the whole region of the recirculation zone, and are presented as pairs of positive and negative components. Reconstructing the 3-D flow field by superimposing the most unstable mode on the base flow shows that the primary vortex is corrugated and the reattachment line exhibits a zigzag shape. Streamwise vortices occur in the primary vortex, accompanied by the formation of a spanwise array of pairs of positive and negative streamwise perturbations. Furthermore, the energy budget analysis suggests that the instability of supersonic BFS flow is probably of the elliptic type.

A step height of h = 7 mm is selected to conduct a DNS study. The spanwise length of the computational region for DNS is four times the wavelength of the most unstable mode. The exponential linear growth rate obtained from DNS shows good agreement with the predicted rate obtained from GSA, and the most unstable mode (mode 1 predicted by GSA) dominates the linear growth stage, in which the contours of the spanwise velocity coincide with the corresponding spanwise velocity perturbations from GSA. Mode 2 emerges in the later linear growth stage and then affects the flow structure together with mode 1. The saturated flow exhibits a broadband low-frequency unsteadiness characteristic. The boundary-layer streaks and surface friction streaks appear in an oscillatory form with a low frequency. In addition, the separation zone is severely distorted, especially the secondary separation zone that fragments into several small bubbles. No transition is observed in the present case.

In conclusion, the numerical simulations and GSA illustrate that supersonic BFS flow is intrinsically unstable to 3-D perturbations when the step exceeds a certain value. The streak-like flow structures that were observed in supersonic BFS flow can be attributed to the global instability of the separated flow.