On the synchronisation of three-dimensional shock layer and laminar separation bubble instabilities in hypersonic flow over a double wedge

Abstract

Linear three-dimensional instability is studied in the shock layer and the laminar separation bubble (LSB) induced by shock-wave/boundary-layer interactions in a Mach 7 flow of nitrogen over a double wedge with a $30^{\circ }\text {--}55^{\circ }$ cross-sectional profile. At a free-stream unit Reynolds number $Re=5.2\times 10^{4}\,{\rm m}^{-1}$ this flow exhibits rarefaction effects and has shock thicknesses comparable to the thickness of the boundary layer at separation. Flow features have been fully resolved using a high-fidelity massively parallel implementation of the direct simulation Monte Carlo method that captures the flow evolution from the inception of three-dimensionality, through linear growth of instabilities, to the early stages of nonlinear saturation. It is shown that the LSB sustains self-excited, small-amplitude perturbations that originate past the primary separation line and lead to spanwise-periodic wall striations inside the bubble and downstream of the primary reattachment line, as known from earlier experiments, simulations and instability analyses. A spanwise-periodic instability, synchronised with that in the separation zone, is identified herein for the first time, which exists in the internal structure of the separation and detached shock layers, and manifests itself as spanwise-periodic cats-eyes patterns in the global mode amplitude functions. The growth rate and the spanwise-periodicity length of linear disturbances in the shock layers and the LSB are found to be identical. Linear amplification of the most unstable three-dimensional flow perturbations leads to synchronised low-frequency unsteadiness of the triple point, with a Strouhal number of $St\approx 0.028$.

1. Introduction

Linear instability in the laminar separation bubble generated by shock-wave/laminar boundary-layer interaction at supersonic and hypersonic speeds has been the subject of intense experimental and, more recently, theoretical/numerical investigations, owing to the ubiquitous nature of this phenomenon and its potential impact on the design of high-speed vehicles. The appearance of regularly spaced streamwise-aligned vortices, or striations, in the reattachment region of nominally two-dimensional supersonic and hypersonic flow over a planar compression corner has been first reported over half a century ago by Ginoux (1958), who went on to relate the vortex spacing with the attached boundary-layer thickness ahead of the corner and document the quantitative effect that the Mach number has on the appearance and spacing of the vortices (Ginoux 1965a,b, 1969, 1971). To date, the interaction of the three-dimensional shock-wave system associated with the compression corner, also first reported experimentally at around the same time by Holden (1963), has not been associated with the process of formation of the striations and motivates the present contribution. However, before turning our attention to the description of oscillations and spanwise modulation in the shock layers, and their association with three-dimensional global instability in the laminar separation bubble (LSB), a brief exposition of our current understanding of the physical origins of the streamwise vortices in the compression corner is warranted, in order to set the scene for the present analysis.

Striations in a compression corner generated a long-standing debate regarding the origins of this phenomenon being attributed to (imperfections at) the leading edge of the flat plate preceding the corner (Simeonides 1992; Simeonides & Haase 1995) and entering the separation region (Navarro-Martinez & Tutty 2005) or an intrinsic flow instability mechanism associated with self-excitation of the LSB formed at the corner itself. In recent years, evidence has been amassed in support of the latter mechanism: fully resolved three-dimensional direct numerical simulations by Shvedchenko (2009) and Egorov, Neiland & Shvedchenko (2011) over a wide range of Mach numbers demonstrated exponential growth of small-amplitude perturbations, the latter taking the form of spanwise-periodic modulation of the wall shear and wall heat transfer footprints of the flow; such striations were found to originate inside the LSB formed at the compression corner and extend beyond the primary reattachment line on the ramp wall. Amplification of the spanwise-periodic perturbations was documented beyond a certain value of the corner angle, scaled according to triple-deck theory (Neiland 1969; Stewartson & Williams 1969; Rizzetta, Burggraf & Jenson 1978a; Neiland et al. 2008). More recently, global linear instability analyses have unequivocally demonstrated the existence and linear amplification of stationary three-dimensional self-excited perturbations originating inside LSBs formed on account of shock impingement on a laminar boundary layer (Robinet 2007; Nichols et al. 2017) or that forming in the steady laminar nominally two-dimensional compression ramp flow (Dwivedi et al. 2019; Cao et al. 2021; Hao et al. 2021).

In the authors’ view, the debate regarding the origins of striations need not lead to the search for mutually exclusive instability paths. Fundamental studies of linear global instability of steady LSBs in the incompressible limit have demonstrated that bubbles can act as amplifiers which promote exponential growth of incoming disturbances, but can also become self-excited due to exponential amplification of intrinsic stationary three-dimensional global modes through an oscillator instability mechanism (Chomaz, Huerre & Redekopp 1988; Theofilis, Hein & Dallmann 2000; Chomaz 2005). Compressibility may alter the relative significance of the two scenarios quantitatively, but not qualitatively. In this context, the interpretation of the experiments of Simeonides (1992) and the direct numerical simulations (DNS) work of Navarro-Martinez & Tutty (2005) would fall in the first category of a bubble acting as an amplifier, while the DNS of Shvedchenko (2009) and Egorov et al. (2011), the global instability analysis of Dwivedi et al. (2019) as well as the combined analysis and experiments of Dwivedi et al. (2020) and Hao et al. (2021) on the compression corner and the global instability analyses of Robinet (2007), Nichols et al. (2017) and Hildebrand et al. (2018) in the related problem of shock-generated laminar separation on a flat plate would be examples of the oscillator scenario.

The present contribution is concerned with the relatively unexplored question of the dynamic behaviour of the shock layer which drives the interaction with the boundary layer and the subsequent appearance of striations downstream of the primary separation line. Shock unsteadiness and the resulting small-amplitude harmonic spatial corrugations on the shock, generated as a result of flow instability, were first discussed in isolation from the boundary layer by Carrier (1949), who addressed the inviscid flow limit and employed linearised Euler equations and Rankine–Hugoniot conditions to study shock interaction with acoustic perturbations. Moore (1954) and Ribner (1954) also monitored an isolated straight shock wave and used inviscid equations to analyse pressure and vorticity waves interacting with the shock, while McKenzie & Westphal (1968) quantified, in terms of the free-stream flow Mach number and the shock angle, the amplitude of acoustic perturbations emitted as a consequence of pressure, vorticity or entropy perturbations traversing a plane stationary isolated shock. Duck & Balakumar (1992) introduced viscosity in their study of self-excitation of a finite-thickness steady normal shock, the latter computed by imposing constant upstream and Rankine–Hugoniot downstream conditions in the framework of the base flow computation methodology proposed by Gilbarg & Paolucci (1953). Their modal analysis revealed that the eigenspectrum only contains (several) branches of damped continuous modes, while their solution of the initial-value problem did not yield reliable large-time response of the shock to incoming perturbations, due to the relatively rich structure of the continuous spectrum. Chang, Malik & Hussaini (1990), Esfahanian (1991) and Malik & Anderson (1991) solved the linear stability eigenvalue problem in the compressible (attached) flat-plate boundary layer (Mack 1984) including far-field boundary conditions based on the Rankine–Hugoniot equations, to analyse the effect of a shock that is locally parallel at a finite distance to the flat plate. Stuckert & Reed (1994) solved the same problem on a cone at very high Mach number, including chemistry effects. Triple-deck theory (Stewartson & Williams 1969; Smith 1986; Neiland et al. 2008) has been employed to understand separated laminar boundary-layer instability in supersonic and hypersonic flows over compression ramps at moderate to high Reynolds numbers (Rizzetta, Burggraf & Jenson 1978b; Cowley & Hall 1990; Smith & Khorrami 1991; Cassel, Ruban & Walker 1995; Korolev, Gajjar & Ruban 2002; Fletcher, Ruban & Walker 2004), however, instability of the shock layer, or its coupling with instability in the boundary layer, are beyond the scope of this theory.

On the other hand, it has long been known (Liepmann, Narasimha & Chahine 1962, 1966) that predictions of the shock structure based on the Navier–Stokes equations increasingly deviate from those delivered by kinetic theory as the degree of rarefaction increases in critical zones such as the internal structure of the shock layer and the Knudsen layer in the vicinity of a solid surface. Liepmann et al. (1962) computed shock profiles using the Bhatnagar–Gross–Krook model of the Boltzmann equation and documented the systematic departure from Navier–Stokes predictions in the low-pressure region of the shock layer, i.e. up to the location of maximum stress, as the free-stream Mach number increases. It then follows that, if the internal structure of the shock is to be resolved, rather than modelled, a numerical approach based on kinetic theory must be used.

Such an effort has been reported by Tumuklu, Levin & Theofilis (2018a); Tumuklu, Theofilis & Levin (2018b), who performed direct simulation Monte Carlo (DSMC) simulations of Mach 16 nitrogen flow over an axisymmetric double cone configuration at unit Reynolds numbers $Re=0.935\text {--}3.74\times 10^{5}\,{\rm m}^{-1}$ and, within the limitations of a two-dimensional simulation, demonstrated a strong coupling between oscillations of the triple point and instability of the laminar separated flow region. The amplitude functions of the least-damped global mode computed in these works comprised structures inside the LSB as well as at the shock system, and included the multiple reflections leading to the lambda-shocklet system between the wall and the slip line downstream of the reattachment point; these features of the global mode became increasingly evident as the Reynolds number at which simulations were performed was increased. In the axisymmetric limit addressed, a steady state was reached after this least-damped global mode decayed linearly, at a constant rate quantified by processing the DSMC signal. Sawant, Levin & Theofilis (2021a,b) studied kinetic fluctuations in the internal non-equilibrium region of one-dimensional shock layers resolved using DSMC and showed that the dominant frequencies are characterised by a Strouhal number, St, defined based on the shock thickness that is two orders of magnitude lower in the shock ($St_{{shock}}\sim O(0.01)$) than that in the upstream ($St_{{shock}}\sim O(1)$). More recently, Klothakis et al. (2022) carried out DSMC simulations in which the temporal and spatial development of harmonic perturbations introduced in a hypersonic flat-plate boundary layer were quantified, before applying linear stability analysis to document the synchronisation of oscillatory perturbations propagating along the leading-edge shock layer with damped perturbations introduced through a wall orifice into the boundary layer.

The present work addresses the three-dimensional analogue of the configuration studied by Tumuklu, Levin & Theofilis (2019), in which the cross-section of the double wedge is extruded along the third spatial direction. Three-dimensional DSMC simulations are performed to study linear instability of the interaction between the shock and the three-dimensional LSB formed in a Mach 7 nitrogen flow over a $30^{\circ }$–$55^{\circ }$ double-wedge configuration. Flow conditions are oriented to those of the experiments performed by Swantek & Austin (2015) and Knisely & Austin (2016) at 8 MJ enthalpy. However, in the present simulations, the number density is kept eight times lower than that in the experiments (here $Re=5.2\times 10^{4}\,{\rm m}^{-1}$), to ensure that all DSMC numerical requirements are met. The kinetic simulation results are analysed to extract characteristics of the small-amplitude three-dimensional perturbations that emerge from numerical noise and are organised into coherent structures both inside the LSB as well as in the separation and detached shock layers. Using resources at the limit of present-day massively parallel computing capabilities it has been possible to fully resolve the (comparable in thickness) internal structure of the shock layer and the boundary layer in the separation region and document, for the first time, the genesis, spatial structure and temporal evolution of small-amplitude shock-layer perturbations and their relation to global instability in the laminar separation zone.

The paper is organised as follows: § 2.1 describes the DSMC algorithm, the specific numerical approach, required numerical parameters and collision models. Details on the initialisation of the three-dimensional simulations are provided in § 2.2, while the features of the two-dimensional base flow, recomputed here for consistency, as well as surface rarefaction effects, are discussed in § 2.3. Results are presented in § 3, starting with quantification of the extraction of linear analysis quantities from the DSMC signal in § 3.1. In § 3.2 the presence and synchronous amplification of a stationary three-dimensional global mode in the LSB and the shock layer is documented, with emphasis placed in § 3.3 on the discussion of spanwise-periodic sinusoidal structures arising from global instability inside the internal structure of the shock layer. Section 3.4 reports the emergence of low-frequency unsteadiness in the shock/LSB system as a consequence of the nonlinear saturation of the amplified global mode. A short discussion summarising the main findings is offered in § 4.

2. Flow modelling methodology

2.1. The DSMC algorithm

The particle-based DSMC method (Bird 1994) is employed to fully resolve the three-dimensional unsteady hypersonic flow field in question. DSMC is a computationally powerful tool, based on the decoupling of molecular motion and intermolecular collisions, in which each simulated particle represents a finite number, $F_{n}$, of actual gas molecules. When all numerical parameters are satisfied, the time-accurate flow field represents a solution of the Boltzmann equation of transport of the molecular velocity distribution function, $f(t,\boldsymbol {r},\boldsymbol {v})$ with respect to time $t$ and position vector $\boldsymbol {r}$, as

(2.1)\begin{equation} \frac{\partial f}{\partial t} + (\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}) f + \left( \frac{\boldsymbol{F}}{m} \boldsymbol{\cdot} \boldsymbol{\nabla}_v \right)f = \left[\frac{{\rm d}f}{{\rm d}t}\right]_{coll} \end{equation}

where $\boldsymbol {v}$ is the instantaneous velocity vector, m is the mass, and $\boldsymbol {\nabla }$ and $\boldsymbol {\nabla }_v$ are gradient operators in physical and velocity spaces, respectively. The first, second and third terms on the left-hand side describe the change of $f$ with time, convection of molecules in space, and convection in velocity space as a result of the external conservative forces per unit mass $\boldsymbol {F}/m$, such as gravity or electric field, which are ignored in this work, respectively. The right-hand side term accounts for changes in $f$ in an element of space–velocity phase space due to intermolecular collisions. Vincenti & Kruger (1965) provide a thorough description of the Boltzmann equation.

The DSMC algorithm typically comprises four major steps of particle mapping, selection for collisions, evaluation of collision outcomes and movement. Based on the choice of boundary conditions, particles are introduced, removed or reflected from the domain boundaries and interacted with the embedded surfaces using gas-surface collision models during a timestep. They are then mapped to a mesh that encompasses the flow domain with cell size, $\Delta x$, that is smaller than the local mean free path of molecules, $\lambda$. Particle pairs are selected for collisions based on the appropriate elastic or inelastic collision cross-section and are then assigned post-collisional instantaneous velocities and internal energies in each computational cell. Macroscopic flow quantities of interest such as pressure, velocities and temperature can be derived from the microscopic properties of simulated particles using the statistical relations of kinetic theory. Finally, based on the post-collisional outcome, particles are moved with a timestep, $\Delta t$, that is lower than the local mean collision time, $\tau$. Conveniently, this physics based approach provides accurate modelling of the internal structure of shocks, their mutual interactions and surface rarefaction effects. This warrants the use of DSMC for detailed modelling of shock-wave/boundary-layer interactions (SBLIs), compared with ad hoc techniques of modelling shocks in numerical solutions of the Navier–Stokes equations that fall short of accurately capturing the internal structure of the shock (Liepmann et al. 1962).

However, to obtain a physically meaningful solution of an unsteady SBLI flow, a number of numerical criteria must be satisfied (see Tumuklu et al. 2018a). To overcome this challenge, an octree-based, three-dimensional DSMC solver known as Scalable Unstructured Gas-dynamic Adaptive mesh-Refinement (SUGAR-3D) (see Sawant et al. 2018) has been utilised in the present work. In summary, the code takes advantage of adaptive mesh refinement (AMR) of coarser Cartesian octree cells to achieve spatial fidelity at regions of strong gradients, cutcell algorithms to correctly model gas-surface collisions in the vicinity of embedded surfaces, domain decomposition strategy based on Morton-Z space filling curves and inclusion of thermal non-equilibrium models, all within a message-passing-interface (MPI) environment that enables massively parallel communication. In the octree-based AMR framework, the collision or $C$-mesh is formed from a user-defined, uniform Cartesian grid with cells, known as ‘roots’, that are recursively subdivided into eight parts until the local cell size is smaller than the local mean free path. Note that a subdivision based on the above criterion is performed only if there are at least 32 particles in a collision cell.

The satisfaction of the numerical criteria in the flow over the three-dimensional double wedge depends on the free-stream conditions, given in table 1. Table 2 provides a summary of the key numerical parameters for these simulations which were found to require $\sim$60 billion computational particles and $\sim$4.5 billion collision cells of the adaptively refined $C$-mesh octree grid. Here, the Knudsen number is based on the length of the lower wedge, 50.8 mm, $F_n$ denotes the number of actual gas molecules represented by a computational particle, while the wall surface is fully accommodated (Bird 1994), i.e. isothermal. The appendix of Sawant et al. (2022) describes the DSMC numerical criteria and convergence study for flows presented here. Figure 6(b,c) of Sawant et al. (2022) show that the ratio $\tau /\Delta t$ is greater than unity in the entire flow domain and the smallest collision cells contain at least eight particles, respectively. With respect to $\lambda /\Delta x$, figure 6(a) of Sawant et al. (2022) shows that the ratio is greater than unity in the entire flow domain, including the interior zones of shock structures and regions of triple, separation and reattachment points. This ratio is between 0.9 and 1, however, in close proximity to the corner of two wedges, also known as the hinge, inside the recirculation region, where the expected errors in transport coefficients of shear viscosity and thermal conductivity are within 6.2 and 3.7 %, respectively (see Alexander, Garcia & Alder 1998). Finally, figure 6(e) of Sawant et al. (2022) also shows that the effect of doubling the number of particles produced the same transient behaviour, where the instantaneous, spanwise-averaged, streamwise velocity differed by at most 7 % in the small region where $\lambda /\Delta x=0.9$ and by 1.5 % outside of it.

Table 1. Physical gas parameters.

Table 2. Numerical simulation parameters.

In terms of collision models, the majorant frequency scheme (MFS) of Ivanov & Rogasinsky (1988) derived using the Kac stochastic model for the selection of collision pairs and the variable hard sphere model for elastic collisions are used. Appendix A describes an improved strategy for the computations of maximum collision cross-section used in the MFS scheme for accurate spectral analysis of unsteady flows simulated on adaptively refined grids. For rotational relaxation, the Borgnakke & Larsen (1975) model is employed with rates of Parker (1959) and DSMC correction factors (see Lumpkin, Haas & Boyd 1991; Gimelshein, Gimelshein & Levin 2002) that account for the temperature dependence of the rotational probability. For vibrational relaxation, the semi-empirical expression of Millikan & White (1963) is used with the high-temperature correction of Park (1984).

2.2. Initialisation of three-dimensional flow

The three-dimensional flow studied in this work is initialised using the two-dimensional steady-state solution over a double wedge simulated by Tumuklu et al. (2019). The latter flow is recalculated here and then replicated on every octree cell along the spanwise direction ($Y$) to form the three-dimensional field at the beginning of the unsteady simulation. Note that the Cartesian coordinates of streamwise, spanwise, streamwise-normal directions are designated as $X$, $Y$ and $Z$, respectively. From the inlet boundary at $X=0$, an inward-directed local Maxwellian flow is introduced at an average number density, bulk velocity and temperature of $n_1$, $u_{x,1}$ and $T_{tr,1}$, respectively. For this work, the SUGAR solver is employed with spanwise-periodic boundary condition in the $Y$-direction, as described in Appendix A. The spanwise extent of the current simulation is $L_y=28.8$ mm, which will be shown to be sufficiently wide to capture linearly growing spanwise-periodic flow structures. In subsequent sections, contours and isocontours showing the detailed spanwise-periodic structures are presented with two periodic wavelengths for clarity.

In terms of the computational strategy used, it takes $\sim 50T$ (flow-through time, $T$, defined in the next section) until the onset of linear instability. To minimise data storage, we estimate this time at which the spanwise-periodic flow structures start to emerge such that they can be easily detected in the background of statistical fluctuations in DSMC. We use an equilibrium criterion along the lines of Hadjiconstantinou et al. (2003) to quantify the amount of spanwise fluctuations about the two-dimensional base state. For example, the standard deviation in the bulk velocity at equilibrium can be calculated as, $\sigma _u=\sqrt {R\langle T_{tr}\rangle _s/\langle N\rangle _s}$, where $\langle T_{tr}\rangle _s$ and $\langle N\rangle _s$ are local spanwise-averaged translational temperature and number of particles, respectively, and $R$ is the gas constant. The linear instability can be easily detected, when the DSMC-computed standard deviation is larger than the above criterion.

Finally, to accommodate billions of computational particles and collision cells, 19.2 k MPI processors were used in this simulation. The spanwise-periodic simulation takes $\sim$5.5 h per through-flow time ($T$) using Intel Xeon Platinum 8280 (‘Cascade Lake’) processors of the Frontera supercomputer (2019) and $\sim$14 h per flow time using Intel Xeon Platinum 8160 (‘Skylake’) processors of the Stampede2 supercomputer (2019). The overall cost of the simulation up to $T=190$ discussed in what follows required $\sim$870 k node hours ($\sim$43.8 million core hours) of computing time.

2.3. Main features of the two-dimensional base flow

Starting with the two-dimensional base flow, figure 1(a) shows the complex features of an Edney-IV type SBLI system (Edney 1968; Babinsky & Harvey 2011). The base flow results are obtained by spanwise and temporally averaging the DSMC flow fields between $T=48$ and 60, where $T$ is the non-dimensional flow time, defined as the time it takes for the flow to traverse the distance $L_s=40$ mm at a free-stream velocity of $u_{x,1}$; $L_s$ is defined as the straight-line distance from the separation point, $P_s$, to the reattachment point, $P_R$, which mark the origin and end of the shear layer inside the separation bubble, respectively. The SBLI features include the interactions of the separation shock with the leading-edge oblique and detached bow shocks at triple points $T_1$ and $T_2$, respectively. Two contact surfaces, $C_1$ and $C_2$, are formed downstream of triple points $T_1$ and $T_2$, respectively, where $C_1$ is between two supersonic streams downstream of the separation shock, and $C_2$ is between the supersonic and subsonic streams downstream of the separation and detached shocks, respectively. For details about the earlier time evolution of the two-dimensional SBLI features over the double wedge, see Tumuklu et al. (2019).

Figure 1. (a) Magnitude of mass density gradient of the base flow, $|\boldsymbol {\nabla }\rho _b|$, normalised by $\rho _1L_s^{-1}$, where $\rho _1=n_1 m$ is free-stream mass density. (b) Overlaid on the image shown in (a) are the wall-normal directions $U$, $S$ and $R$, and the numerical probes $a, b, c, d, r, s$ and $t$, defined in the text.

The image shown in figure 1(a) is repeated in figure 1(b), in order to define a number of probes that will aid the discussion of the next sections. Probes are placed at locations $a$: at the foot of the separation shock, ($X=32.817$ mm, $Z=13.285$ mm), $b$: inside the LSB, ($X=48.496$ mm, $Z=24.270$ mm), $c$: at the contact surface, ($X=65.191$ mm, $Z=56.593$ mm), $d$: on the detached shock, ($X=52.860$ mm, $Z=59.226$ mm), $r$: at reattachment, ($X=64.396$ mm, $Z=44.358$ mm), $s$: on the separation shock, ($X=44.165$ mm, $Z=32.597$ mm) and $t$: at the triple point $T_2$, ($X=48.347$ mm, $Z=41.624$ mm). Probes $a$ and $d$ define lines $a\text {--}t$ and $t\text {--}d$, denoting the intersection with the $OXZ$ plane of two planes passing through the separation and detached shocks, respectively. Wall-normal directions $U$, $S$ and $R$ are respectively defined at the foot of the separation shock and either side of the triple point, $T_2$. The vertical projection of the foot of the arrows $U$, $S$ and $R$ on the X-axis is at 32.5, 49 and 66.1 mm, respectively. From this point onward, the DSMC-derived instantaneous profiles shown are noise filtered using the proper orthogonal decomposition (POD) procedure discussed in Appendix B.

The use of a particle-based approach permits analysis of surface rarefaction effects; results corresponding to the base flow are shown in figure 2(a) for the local streamwise velocity slip, $V_{s}$, normalised by $u_{x,1}$, while the local mean free path adjacent to the wall, $\lambda$, as well as the translational temperature jump at the surface are presented in figure 2(b), respectively normalised by free-stream mean free path $\lambda _1$ and $T_{tr,1}$.

Figure 2. Surface macroscopic flow quantities in the base state, where the profiles are time averaged from $T=48$ to 53. (a) Surface velocity slip in the local streamwise direction, $V_{s}$. (b) Value of $\lambda$ adjacent to the wall and temperature jump $T_{s}$. (c) The heat transfer and pressure coefficients, $C_h$ and $C_p$, respectively.

Figure 2(a) shows a maximum $V_s$ equal to 2.16 % of $u_{x,1}$ at the leading edge ($X=10$ mm), similar to the value of 2.45 % obtained by Tumuklu et al. (2019) in a two-dimensional flow simulation. The large slip at the leading edge is due to the increased rarefaction of gas induced by steep gradients of the leading-edge shock. Here, $V_s$ decreases along the local streamwise direction to 0.6 % at $X=32$ mm, similar to the decrease of $\lambda$ adjacent to the wall as seen from figure 2(b). From $X=32$ to 36 mm in the interaction region of the separation shock with the boundary layer, $V_{s}$ rapidly decreases to zero at the separation point, $P_S$. Inside the recirculation zone, from $P_S$ to $P_R$, the point of reattachment, $V_{s}$ is negative because the flow impinging on the wall is opposite to the local streamwise direction; $V_{s}$ and $\lambda$ remain constant on the lower wedge, where the latter is approximately 3.69 % of the free-stream mean free path, $\lambda _1$. On the upper slant surface, $V_{s}$ and $\lambda$ both increase such that the rate of increase is significantly larger downstream of reattachment while $T_s$ follows a similar variation as $V_{s}$ and $\lambda$, as seen from figure 2(b). Velocity slip and temperature jump are rarefaction effects that are proportional to the size of the Knudsen layer in the vicinity of the wall (Chambre & Schaaf 1961; Kogan 1969). Since the Knudsen layer is approximately of the order of $\lambda$, which is inversely proportional to number density, $n$, and directly proportional to the translational temperature, $T_{tr}^{\omega -0.5}$, where $\omega =0.745$ is the viscosity index of the gas, the profiles of $V_s$ and $T_s$ are similar to the profile of $\lambda$ in the vicinity of the surface, for our assumption of a fully diffuse surface.

The variations in surface heat flux and pressure coefficients, $C_h$ and $C_p$, respectively, are shown in figure 2(c). The change in $C_h$ along the local streamwise direction is similar to $V_s$ and $T_s$ and $C_p$ is constant on the lower wedge, except for a sharp increase in the vicinity of the separation point from $X=32$ to 38 mm. Inside the recirculation zone on the lower wedge, $C_p$ plateaus, but on the upper wedge it increases rapidly up to the corner of the wedge. It must be noted that, although figure 2(c) shows a small dip in $C_p$ at the hinge equal to 1.36 % of the plateau at $X\sim 50$ mm, we do not observe secondary vortices in the separation region. This is consistent with the value of the scaled angle $\alpha \sim 2$ of our simulations, which is much smaller than that required for secondary vortices to emerge ($4 \le \alpha \le 5$) (see Shvedchenko 2009). The scaled angle is calculated based on the definition used by Gai & Khraibut (2019) (eq. (1.1)), which was originally defined by Stewartson (1970) and Rizzetta et al. (1978a).

3. Three-dimensional linear instability mechanisms

3.1. Analysis of the three-dimensional DSMC signal

During the early time of the simulation, when deviations from the imposed steady state are still small in magnitude, perturbations $\tilde {Q}$ of an unsteady macroscopic flow quantity $Q = (n,u_x,u_y, u_z, T_{tr}, T_{rot}, T_{vib})^{\text {T}}$ may be obtained by subtracting from the unsteady three-dimensional full DSMC field the steady, two-dimensional base flow state $\bar {Q} = (\bar {n},\bar {u}_x,0, \bar {u}_z, \bar {T}_{tr}, \bar {T}_{rot}, \bar {T}_{vib})^{\text {T}}$ computed in § 2.2 and imposed as initial condition in the simulation

(3.1)\begin{equation} \tilde{Q}(x,y,z,t)=Q(x,y,z,t)-\bar{Q}(x,z). \end{equation}

In the disturbance field vector $\tilde {Q}= (\tilde {n},\tilde {u}_x,\tilde {u}_y, \tilde {u}_z, \tilde {T}_{tr}, \tilde {T}_{rot}, \tilde {T}_{vib})^{\text {T}}$, $\tilde {n}$ denotes the perturbation number density, $\tilde {u}_x,\tilde {u}_y, \tilde {u}_z$ are perturbation velocities in the $X$, $Y$ and $Z$ directions and $\tilde {T}_{tr}, \tilde {T}_{rot}, \tilde {T}_{vib}$ are perturbation translational, rotational and vibrational temperatures, respectively. The smallness of the perturbations at early times in the simulation permits decomposition of the DSMC-computed perturbation flow fields $\tilde {Q}$ according to linear stability theory into a spanwise-independent vector $\hat {Q}(x,z)$, comprising the two-dimensional amplitude functions, and a phase function, $\varTheta$, according to (Theofilis 2000; Theofilis & Colonius 2003)

(3.2)\begin{equation} \tilde{Q}(x,y,z,t) = \hat{Q}(x,z)\exp{(\text{i}\varTheta)} + \text{c.c.}, \end{equation}

where

(3.3)\begin{equation} \varTheta = \beta y - \varOmega t, \end{equation}

Here, $\beta =2{\rm \pi} /L_y$ is the real spatial wavenumber corresponding to the spanwise wavelength, $L_y$, of the global mode, ${\varOmega =\varOmega _r + \text {i} \varOmega _i}$ is a complex quantity, whose real part, $\varOmega _r$, indicates the frequency and the imaginary part, $\varOmega _i$, is the temporal growth rate and $\text {c.c.}$ indicates complex conjugation so that $\tilde {Q}$ is real. The sign convention in (3.3) indicates that $\varOmega _i>0$ signifies exponentially growing linear perturbations.

3.2. Linear instability in the three-dimensional shock layer/LSB interaction region

Analysis of the DSMC results was performed during two time windows, $30\le T\le 68$ and $68 \le T\le 190$. The lower limit of the first bracket was chosen after initial transients in the solution had subsided, while the upper limit in the second time interval was chosen in order to study the long-time development of the instability, once nonlinear saturation had been reached. A qualitative change that occurs after $T=68$ in the shock structure and its long-time ($T>110$) effect on the LSB will be discussed in § 3.4; in this and the next section an in-depth discussion of the results obtained during the linear growth of perturbations is presented first.

The temporal evolution of the spanwise perturbation velocity component, $\tilde {u}_y$, normalised by $u_{x,1}$, within the time of linear growth, indicated by the greyed out box in figure 3, is presented at probe $s$ inside the shock and probe $b$ inside the LSB, both at two spanwise locations, $C$ and $D$, to be defined shortly. The immediate observation made is that, both inside the shock layer and in the LSB, $\tilde {u}_y$ grows exponentially and, up to this time, monotonically. In order to verify the existence of an underlying linear instability mechanism and quantify its parameters, a two-dimensional linear fit of the full three-dimensional field of the perturbations is performed according to (3.2) using the generalised least-squares method of the Python LMFIT (Version 1.0.1) function. The results returned for the mean value of the unknown fit parameters, $\varOmega _i$, $\hat {Q}$, $\varOmega _r$, alongside the $1\sigma$-uncertainty (standard error) of these estimates, are presented in table 3, where the growth rates obtained for macroscopic flow quantities at probes $b, s, r$ and $c$ defined in figure 1 are shown.

Figure 3. Temporal evolution of spanwise perturbation velocity, (a) in the separation shock (probe $s$) and (b) the bubble (probe $b$), at spanwise locations $C$ ($y/L_y=0.13$) and $D$ ($y/L_y=0.63$).

Table 3. Two-dimensional linear curve fit parameters obtained using (3.2) and (3.3). Fit parameters are obtained at probe $b$ unless explicitly indicated otherwise.

The results of curve fits at all locations examined indicate the presence of an amplified stationary global linear instability. The average growth rate of each flow quantity listed in the table is $\varOmega _i=5.0$ kHz, with a standard deviation of $\sim$ 6.7 %. Significantly, the absolute magnitude of the small amplitude of spanwise perturbation velocity, $\tilde {u}_y$, is of the same order at probes $r$, $b$, $s$, and an order of magnitude lower at probe $c$: the maximum deviation of $\sim$11.4 % is found at probe $c$, which is outside of the region of shock and LSB interaction.

The temporal evolution of $\tilde {u}_y$ at $50 \le T \le 90$ over the entire span of the separation bubble is shown in figure 4(a) as raw DSMC data, while figure 4(b) shows a two-dimensional linear fit of the data shown in figure 4(a) using (3.2) and (3.3). A well-defined harmonic pattern of $\tilde {u}_y$ emerges along the span and the choice of points $C$ and $D$ becomes clear, as they represent the spanwise locations at which $\tilde {u}_y$ attains a local maximum and minimum, respectively.

Figure 4. Temporal evolution of spanwise perturbation velocity, $\tilde {u}_y$, inside the LSB, along the entire span. (a) Raw DSMC data. (b) Two-dimensional linear fit of the result shown in (a).

The emergence of linear instability in the DSMC simulation is demonstrated in figure 5, where the probe results shown in figure 3 at location C are plotted again in semi-logarithmic scale; superposed upon the DSMC results a straight line is plotted, having a slope $\varOmega _i=5.0$ kHz, as indicated by the average value of 5 kHz of the amplification rate calculated earlier. A number of significant conclusions can be drawn on the basis of these results. First, perturbations grow exponentially by at least one order of magnitude inside the separation bubble and the shock layer, up to times $T\approx 110$ and $T\approx 130$ at probes $b$ and $s$, respectively. Second, within the one standard deviation error bar of 0.335 kHz (6.7 %) calculated earlier, growth of perturbations occurs at the same amplification rate inside the bubble and the shock, as indicated by the parallel slopes of the straight dashed line and those of the probe data. Third, at these conditions nonlinear saturation, indicated by the oscillatory evolution of the signal that will be discussed in the next section, occurs first inside the LSB and sets in later in the shock layer.

Figure 5. Same data at location $C$ as in figure 3, plotted in a semi-logarithmic scale. The dashed line shows the average value of the growth rate shown in table 3, $\varOmega _i\approx 5$ kHz, i.e. $\varOmega _i L_s/u_{x,1}\approx 5.247\times 10^{-2}$. The colour for the data at probe $b$ has been changed to distinguish from the data at probe $s$.

To the best of the authors’ knowledge, this is the first time that (unsteady) DSMC simulations have captured the growth of linear instability in the context of three-dimensional DSMC simulations. The essential qualitative difference between the present three-dimensional and the earlier two-dimensional results of Tumuklu et al. (2019) is the discovery of unstable three-dimensional perturbations in the three-dimensional configuration, as opposed to the damped eigenmodes in the two-dimensional work, that ultimately led to a steady state being reached in the two-dimensional flow; the present results show that no such steady state exists in the three-dimensional double-wedge flow at these parameters. Notwithstanding differences in Mach number, a comparison of the present results with those of instability analysis in the double wedge by Sidharth et al. (2018), who analysed a Mach 5 hypersonic flow of a calorically perfect gas, reveals that the present non-dimensional, average growth rate $\varOmega _i=0.0057$ (growth rate results are non-dimensionalised by multiplying $\varOmega _i$ by the boundary-layer thickness at separation, $\delta _{99}=3.35$ mm, and dividing by the free-stream velocity downstream of the leading-edge shock, $u_{x,2}=2930.8\,{\rm m}\,{\rm s}^{-1}$, the latter obtained from inviscid shock theory (Anderson 2003) for the observed shock angle of 41$^{\circ }$) is nearly an order of magnitude larger than the $\varOmega _i=7.5\times 10^{-4}$ obtained in the referenced work. This can be explained by the difference in the angles between the upper and lower wedges, which is $\Delta \theta =15^{\circ }$ in the present and $\Delta \theta =8^{\circ }$ in the referenced work and is consistent with the strong destabilising effect of $\Delta \theta$ found by Sidharth et al. (2018), and also known from the related compression ramp flow (e.g. Cao et al. 2021). Finally, linear instability was found not to affect surface flow quantities. Spanwise variations of the streamwise velocity slip and the temperature jump at the wall were found to be negligible during linear growth of disturbances, such that these quantities are the same as the base state profiles shown in figures 2(a) and 2(b), respectively. The spanwise velocity slip (not shown) was also found to be negligibly small with a maximum spanwise variation of only 0.078 % of $u_{x,1}$.

The question of the spatial origin of linear instability and the spanwise periodicity seen in figure 4 is addressed next. Two locations are selected, $U$ and $S$, as shown in figure 1(b), located upstream of the separation shock and in the separation region, respectively. Figure 6(a) shows the variation of the local streamwise velocity, $u_{t,l}$ along the local wall-normal direction, subscript ‘$t$’ standing for the local streamwise (wall-tangential) component of velocity and ‘$l$’ referring to the lower wedge surface, as a function the wall-normal height $H_l$. Figure 6(b) shows the absolute maximum spanwise variation of $u_{t,l}$ during linear perturbation growth, calculated based on the difference in $u_{t,l}$ at two locations, $A$, at $Y/L_y=0.88$ and $B$, at $Y/L_y=1.38$ of the spanwise sinusoidal modulation induced by linear instability; the significance of these locations $A$ and $B$ will be further elucidated in the next section.

Figure 6. (a) Boundary-layer profiles in the base state along wall-normal directions $U$ and $S$. Local streamwise velocity, $u_{t,l}$ is normalised by $u_{x,1}$, $H_l$ being the distance along the wall-normal direction. (b) The absolute maximum spanwise variation in $u_{t,l}$ along $B$ and $S$.

Common characteristics of the boundary-layer profiles at $U$ and $S$ are the respective non-zero streamwise velocity component at the wall, due to surface rarefaction effects discussed in § 2.3, and the generalised inflection points (GIPs), indicated by open circles. The essential difference between the profiles at the two locations concerns the respective maximum spanwise variation: at $U$, this is only within 0.05 % of $u_{x,1}$, indicating that the flow upstream of the separation shock is essentially two-dimensional. By contrast, at $S$ the maximum spanwise variation peaks at the inflection point, $H_l/L_y=0.113$, where it is equal to 1.8 % of $u_{x,1}$. In addition, a second local maximum of 0.72 % of $u_{x,1}$ is observed inside the separation shock layer at $H_l/L_y=0.40$. Attention is thus turned to correlating this spanwise inhomogeneity inside the shock layer with global instability in the LSB.

3.3. Correlation of linear instability in the shock layer and the LSB

Figures 7(a) and 7(b) present probe data of perturbation quantities along specific wall-normal and spanwise field lines, respectively. It should be stressed here that the flow components shown are representative of all perturbation quantities, all of which exhibit the same qualitative behaviour, but are not shown for brevity. Figure 7(a) shows the variation of the gradient magnitude of the pressure perturbation, $|\boldsymbol {\nabla } \tilde {p}|$, normalised by $p_1L_s^{-1}$ ($p_1$ being the free-stream pressure) as a function of wall-normal distance, $H_l$, along the $S$-direction. Data are shown at the same two spanwise locations, $A$ and $B$, the choice of which will be discussed shortly. The rapid increase of $|\boldsymbol {\nabla } \tilde {p}|$ at $H_l=0.36L_y$ is indicative of the separation shock, inside of which the value of $|\boldsymbol {\nabla } \tilde {p}|$ far exceeds that in the vicinity of the surface. This criterion was used to indicate the separation shock region in figure 6(a). It is worth pointing out here that the measured thickness of the shock layer based on this criterion, $0.083L_y=2.39$ mm, is comparable to the boundary-layer thickness at separation, $\delta _{99}=3.35$ mm.

Figure 7. (a) Wall-normal probe data for $|\boldsymbol {\nabla } \tilde {p}|$, along direction $S$ as a function of wall-normal height at two spanwise locations $A$ and $B$ defined in the legend. (b) Spanwise probe data for $\tilde {u}_{t,l}$ inside the separation shock at $H_l/L_y=0.4$ on plane $S$.

The choice of spanwise locations $A$ and $B$ can be understood from results shown in figure 7(b). Here, the spanwise variation of $\tilde {u}_{t,l}$ inside the separation shock layer is plotted at a wall-normal height $H_l/L_y=0.4$ along the $S$-direction. It can be clearly seen that the spanwise locations $A$ and $B$ correspond to the trough and peak of a sinusoidal modulation of $\tilde {u}_{t,l}$.

The spanwise structure of small-amplitude perturbations inside the separation and detached shock layers, alongside that inside the LSB, is shown in figures 8(a) and 8(c), respectively corresponding to wall-normal planes at the locations $S$ and $R$ defined in figure 1(b). Figures 8(b) and 8(d) show the same spanwise perturbation velocity component on the separation shock plane $a\text {--}t$ and the detached shock plane $t\text {--}d$ defined in figure 1(b). In all four of these figures contours of the spanwise perturbation velocity component are shown. This quantity is zero at the start of the simulation and only arises on account of linear instability as time progresses. On figures 8(a) and 8(c) the approximate boundaries of the separation and detached shock regions, as well as the envelope of the LSB and the location of the contact surface are marked by dashed/dash-dotted horizontal lines.

Figure 8. Field data of the spanwise perturbation velocity component, $\tilde {u}_y$, normalised by $u_{x,1}$ during linear growth of perturbations, on the planes $S$ in (a), $a\text {--}t$ in (b), $R$ in (c) and $t\text {--}d$ in (d) as denoted in figure 1(b). Overlaid line contours: (black dashed line) $|\boldsymbol {\nabla } \tilde {p}|$ and (black dashed dotted line) $\tilde {\omega }_y=0$.

On the $S$-plane in figure 8(a), spanwise-periodic flow structures are seen inside the separation bubble between the wedge surface, $H_l=0$, and the outer envelope of the separation bubble, a $H_l\approx 0.15L_y$, where the spanwise vorticity, $\tilde {\omega }_y$, is zero. The outer envelope of the bubble can also be seen to be corrugated sinusoidally, on account of the stationary growing LSB instability, the three-dimensional footprint of which leads to the well-known striations that originate inside the LSB (Ginoux 1958; Shvedchenko 2009; Dwivedi et al. 2019; Hao et al. 2021). This spanwise-periodic global instability is identified here for the first time in the context of kinetic theory simulations.

Further away from the wall, figure 8(a) shows, also for the first time, that spanwise-periodic flow structures form inside the strong gradient region of the separation shock layer, at a height $0.36 < H_l/L_y < 0.44$ away from the wall. These structures are in phase with those inside the LSB and have the same spanwise-periodicity length. This result is consistent with the maximum spanwise variation seen in the probe data and the DSMC signal post-processing discussed in the previous section, which revealed nearly identical growth rates inside the LSB (probe $b$) and the separation shock (probe $s$). The conjecture put forward based on interpretation of the probe data can now be further strengthened on the basis of the field data shown, namely that self-excited linear instability leads to spanwise-periodic structures in perturbation flow quantities within the shock layer, with a spanwise wavelength of $L_y$, identical with that seen in the LSB.

On the cut plane $a\text {--}t$ in figure 8(b), it is seen that the spanwise-periodic structures are present in the entire separation shock layer from the foot of the separation to the triple point $T_2$. The region of interaction of the oblique shock with the separation shock is marked on the $a\text {--}t$ plane at the distance of 0.474$L_y$ to 0.613$L_y$. In this region, the magnitude of $\tilde {u}_y$ increases and striation patterns emerge further along the separation shock.

On the $R$-plane, downstream of the triple point, figure 8(c) shows contour plots of the spanwise perturbation velocity component, $\tilde {u}_y$, which also exhibits spanwise-periodic structures inside the reattached boundary layer; together with the result shown in figure 8(a), this result shows that the striations originating inside the LSB extend at least up to the location of plane $R$ on the downstream wedge wall. Along the wall-normal direction, $H_u$, spanwise-periodic structures are also present in the vicinity of contour line $\tilde {\omega }_y=0$ at a height of $H_u=0.36L_y$, in the vicinity of contact surface $C_2$ downstream of the triple point $T_2$, which clearly points to the contact surface being corrugated by perturbations that are in phase with the striations on the wall. Further away from the wall, the contour lines of $|\boldsymbol {\nabla } \tilde {p}|$ at $H_u=0.61L_y$ and $0.677L_y$ indicate the approximate layer of the detached shock. Inside this layer the spanwise-periodic flow structures are still noticeable, but their amplitude has subsided compared with that in the separation shock, indicating that the spanwise-periodic global mode is concentrated in the zone of interaction between the separation and detached shock-layer system and the LSB.

The presence of spanwise-periodic structures in the detached shock layer from the triple point $T_2$ to probe $d$ is seen on the cut plane $t\text {--}d$ in figure 8(d). The peak magnitude of these structures is at least a factor of two lower than the structures in the separation shock and an order of magnitude lower than those in the separation bubble. This shows that the synchronisation of global instability is strongest between the separation shock and the LSB, and the amplitude of the perturbations in the shock layer declines downstream of the triple point. It can further be seen that the structures inside the detached shock layer are 180$^{\circ }$ out of phase with structures in the separation shock, shown earlier on plane $a\text {--}t$.

Returning to the analysis of the spanwise structures inside the shock layer and the separation bubble on the $S$ plane, figure 9 confirms that spanwise-periodic perturbations are present in all perturbation flow quantities, namely number density, $\tilde {n}$, streamwise velocity with respect to the lower wedge (i.e. direction perpendicular to $S$), $u_{t,l}$, wall-normal velocity with respect to the lower wedge (i.e. direction of $S$), $u_{n,l}$, translational temperature, $\tilde {T}_{tr}$, rotational temperature, $\tilde {T}_{rot}$, and vibrational temperature, $\tilde {T}_{vib}$. The minimum (negative) and maximum (positive) values of spanwise structures in $\tilde {u}_{t,l}$, $\tilde {u}_{n,l}$, and $\tilde {n}$ are also at spanwise locations $A$ and $B$, respectively.

Figure 9. Contours of perturbation quantities on the plane defined along $S$. (a) Number density, $\tilde {n}$, (b) streamwise perturbation velocity, $u_{t,l}$, (c) wall-normal velocity, $u_{n,l}$, (d) translational temperature, $\tilde {T}_{tr}$, (e) rotational temperature, $\tilde {T}_{rot}$, and ( f) vibrational temperature, $\tilde {T}_{vib}$, showing the cats-eyes pattern of instability in the shock layer. Number density, velocities and temperatures are normalised by respective upstream values.

Interestingly, closer inspection of the relative magnitude of the perturbation velocity components, all scaled by the same free-stream value, reveals qualitatively different results inside the LSB and at the shock layer. In the LSB, the tangential and spanwise perturbation velocity components have comparable magnitudes, while the wall-normal component is an order of magnitude smaller in relative terms. By contrast, the peak oscillation amplitude in the wall-normal perturbation velocity component is two orders of magnitude higher in the shock layer than inside the LSB ($|\tilde {u}_{n,l}|_{max}\approx 0.3$ inside the shock layer, as opposed to $|\tilde {u}_{n,l}|_{max}\approx 3\times 10^{-3}$ in the bubble), as seen in the inset graphs that magnify all shock-layer results. This behaviour, although not as pronounced, can also be seen in the number density and perturbation temperature components.

Furthermore, all three perturbation temperatures have primary spanwise structures adjacent to the wall having minimum and maximum values at spanwise locations $B$ and $A$, respectively, i.e. they are 180$^{\circ }$ out of phase with perturbation structures of velocities and number density. Perturbations in the normalised $\tilde {T}_{tr}$ and $\tilde {T}_{rot}$ also exhibit secondary structures immediately above the primary structures within $0.1 < H_l < 0.15$. From a qualitative point of view, all perturbation components inside the shock layer are found to feature the same cats-eyes pattern, seen in all insets of figure 9. The degree of thermal non-equilibrium seen in the perturbation temperatures presented in figures 9(d)–9( f) is significant, and originates in differences in the base flow temperatures, shown in figure 10(a), where thermal non-equilibrium is present in the entire region between the shock structure and the wedge surface along the three directions shown. This results in the large differences in the shock, denoted by the blue dashed line at $H_l/L_y = 0.36-0.44$ in figure 10(b), where relatively smaller differences can be seen in the bubble, at $H_l/L_y=0-0.15$, in all three temperatures. Consistent with the faster relaxation between translational and rotational modes, those two temperatures have an analogous spatial dependence.

Figure 10. (a) Base flow temperatures, $\bar {T}_{tr}$, $\bar {T}_{rot}$, $\bar {T}_{vib}$, along $U$, $S$ and $R$. (b) Wall-normal probe data of perturbation temperatures, $\tilde {T}_{tr}$, $\tilde {T}_{rot}$, $\tilde {T}_{vib}$, along $S$ at the spanwise location $A$.

Figure 11 presents three different views of the spanwise perturbation velocity component over the entire calculation domain. In figure 11(a) the perturbation is shown on the $OXZ$ plane and two characteristic cut planes, $P$ and $Q$, defined by their respective normal vectors $[-0.7193\hat {i} + 0.6946\hat {k}]$ and $[0.7193\hat {i} - 0.6946\hat {k}]$. Here, $P$ is taken approximately tangential to the dividing streamline of the LSB, while $Q$ is taken at a distance from the corner, so as to cross the triple point of the steady flow. Figure 11(b) shows three-dimensional isosurfaces of $u_y$ underneath plane $P$, while figure 11(c) presents spanwise-periodic structures in the shear layers downstream of triple points $T_1$ and $T_2$ defined in figure 1(a). All of these results, as well as consistent results seen in the isocontours of vorticity components in Appendix B, all confirm the three-dimensional, spanwise-periodic nature of an instability that encompasses the LSB region in tandem with the separation and detached shock layers and the triple point $T_2$ that connects the latter two.

Figure 11. (a) Isocontours of spanwise perturbation velocity $\tilde {u}_y$ on the $OXZ$ plane and definition of cut planes $P$ and $Q$. (b) Same result for $\tilde {u}_y$, plotted underneath cut plane $P$ and denoting spanwise-periodic striations inside the LSB. (c) Same plotted upstream of cut plane $Q$ and denoting spanwise-periodic modulation in the vicinity of contact surfaces (shear layers) downstream of triple points.

Summarising, it is seen that during linear perturbation growth the steady two-dimensional base flow becomes unstable to a self-excited, small-amplitude, three-dimensional spanwise-periodic stationary global mode. The spanwise perturbation velocity component, $\tilde {u}_y$, which was zero at the beginning of the simulation, attains a sinusoidally varying amplitude not only inside the separation bubble, as known from earlier work (e.g. Ginoux 1958; Theofilis et al. 2000; Shvedchenko 2009; Dwivedi et al. 2020), but also inside the separation and detached shock layers and the shear layers formed downstream of triple points. This result extends the findings of Tumuklu et al. (2018a,b, 2019) in the two-dimensional counterpart of this configuration but, contrary to the two-dimensional limit, in the present three-dimensional environment the observed global instability is unstable with synchronised spanwise-periodic small-amplitude perturbations amplifying exponentially within both the LSB and shock layer. Linear growth of the global mode inside the laminar separation leads to striations on the wall, while in the shock layer a spanwise-periodic cats-eyes pattern is formed in all perturbation components. Ultimately, linear growth will lead to nonlinear saturation and attention is turned next to analysis of the three-dimensional DSMC predictions at later times.

3.4. Low-frequency unsteadiness of the shock and separation bubble

Turning to the dynamics of the linear flow instability at times beyond exponential growth of the stationary spanwise-periodic global mode, attention is paid to low-frequency unsteadiness, that has been reported over a wide range of $Re\sim O(10^{6}-10^{8})$ in a number of numerical solutions of the Navier–Stokes equations (see e.g. Pirozzoli & Grasso 2006; Piponniau et al. 2009; Touber & Sandham 2009; Priebe & Martín 2012; Clemens & Narayanaswamy 2014; Gaitonde 2015), as well as many experiments (Dussauge, Dupont & Debiève 2006), where shock-induced turbulent boundary-layer interactions were studied in a variety of configurations. Low-frequency unsteadiness is identified in this context as that corresponding to Strouhal numbers in the range $0.01 \le St \le 0.05$, while Strouhal numbers within the same range have also been reported for supersonic laminar and transitional shock/boundary-layer interactions (Sansica, Sandham & Hu 2016; Threadgill, Little & Wernz 2019).

A qualitative change in the evolution of the signal commences in the triple point at $T=68$. Figures 12(a) and 12(b) present the spatio-temporal variation of normalised perturbation spanwise velocity, $\tilde {u}_y$, at probe $b$ in the separation bubble, where linearly growing spatial structures start to exhibit sinusoidal oscillations in time with an average period of 46$T$. At probe $t$ near triple point $T_2$, the variation of normalised perturbation number density, $\tilde {n}$, is shown in figures 12(c) and 12(d) where two time scales with periods of 51$T$ and 41$T$ (average period of 46 $T$) can be seen. The corresponding Strouhal number corresponding to these time periods is,

(3.4)\begin{equation} St = \frac{f L_{s}}{u_{x,2}} = 0.0283 \pm 0.003, \end{equation}

which is defined based on the length of the separation bubble, $L_s=40$ mm and the velocity downstream of the leading-edge oblique shock, $u_{x,2}=2930.8$ m s$^{-1}$. This Strouhal number is within the low-frequency range reported in the literature of turbulent SBLI. The triple point starts to oscillate at $T=68$, and its motion remains two-dimensional up to approximately $T=85$, as indicated by the lack of variation in $\tilde {n}$ along the spanwise direction. Past this time, three-dimensional linear instability manifests itself and low-frequency unsteadiness is observed, causing spanwise modulation of structures in $\tilde {n}$ and their oscillations in time.

Figure 12. Low-frequency unsteadiness at probes $b$ and $t$ in the separation bubble and at the triple point, $T_2$, respectively. (a) At probe $b$, the temporal evolution of perturbation spanwise velocity, $\tilde {u}_y$, normalised by $u_{x,1}$. (b) At spanwise locations $C$ ($Y/L_y=0.13$) and $D$ ($Y/L_y=0.63$), the normalised $\tilde {u}_y$ indicating the period of unsteadiness (see the greyed out region). (c) At probe $t$, the temporal evolution of perturbation number density, $\tilde {n}$, normalised by $n_1$. (d) Normalised $\tilde {n}$ at locations $A$ ($Y/L_y=0.88$) and $B$ ($Y/L_y=1.38$) as a function of $T$.

In summary, as time progresses, low-frequency oscillation follows exponential growth of a self-excited stationary three-dimensional global mode in the separation bubble and the shock layer. As far as self-excitation of the laminar separation zone is concerned, the appearance of low-frequency oscillation in the related configuration of laminar separation generated by a shock impinging on a flat-plate geometry has been attributed by Boin et al. (2006) to a supercritical Hopf bifurcation of the flow, following the oscillator scenario discussed by Theofilis et al. (2000) in the incompressible LSB case.

The present kinetic-theory results fully concur with this prediction regarding the existence of a self-excited three-dimensional stationary global instability of the LSB. Moreover, we find that the underlying global mode is not confined to the LSB but includes the shock layer as an integral part. Past linear amplification of this stationary mode, low-frequency synchronised oscillations arise both in the LSB and in the (fully resolved) shock layer. Current computing capabilities have not permitted us to pursue the time integration of the DSMC equations in order to further quantify additional frequency content generated past the initial stages described herein, or follow flow into transition to turbulence; this task will be pursued in future efforts.

4. Summary and conclusions

Linear instability mechanisms in laminar hypersonic flow of nitrogen over a three-dimensional double wedge with a $30^{\circ }\text {--}55^{\circ }$ cross section have been studied using the particle-based DSMC method, as implemented in the massively parallel SUGAR solver (Sawant et al. 2018). The temporal evolution of the DSMC field was computed over a span of $\sim$200 through-flow time units that included the inception of three-dimensionality, linear perturbation growth and nonlinear signal saturation. The fully resolved kinetic approach permitted space- and time-accurate modelling of the internal structure of shocks, surface rarefaction effects, thermal non-equilibrium as well as, for the first time, time-accurate description of small-amplitude three-dimensional self-excited perturbations that arise from interaction of the shock system with the three-dimensional LSB both inside the separation and detached shock layers and in the LSB.

In line with the findings of Tumuklu et al. (2018b) for Mach 16 flow over an axisymmetric double cone, and Tumuklu et al. (2019) for the two-dimensional counterpart of the present double wedge flow, the LSB was found to be strongly coupled with the separation and detached shock system. In contrast to the results in the two-dimensional configuration, linear instability was found to be exponentially amplified in the present three-dimensional context. The origin of instability was traced in the vicinity of primary separation, i.e. the flow region of interaction of the separation shock and the bubble. This finding clearly points out that the observed mechanism is one of self-excitation and not the result of amplification of incoming upstream disturbances, a result which generalises the concept of oscillator, introduced in an incompressible context either in weakly non-parallel (Chomaz et al. 1988; Chomaz 2005) or in fully non-parallel flow (Theofilis et al. 2000), to include the interaction of the entire separation and detached shock system with the three-dimensional LSB.

The presence of amplified linear instability leads to the formation of spanwise-periodic flow structures in three-dimensional perturbations of macroscopic flow quantities inside the LSB, corresponding to the well-known striation wall pattern. More importantly, the present study showed that spanwise-periodic linear instability also exists in the internal structure of the separation and detached shock layers and in the vicinity of contact surfaces downstream of triple points; this shock-layer instability, which takes the form of cats-eyes patterns, is seen for the first time in the present work. Linear instability in the shock layer and the LSB are synchronised, with the spanwise-periodicity length of the structures forming inside the bubble and in the shock system being identical, and their amplitude growing with the same average, linear temporal growth rate, $\varOmega _i=5.0\,{\rm kHz}\pm 0.16\,\%$ at the conditions examined. The slow linear growth and long time scale of low-frequency unsteadiness suggest that experimental test times must be sufficiently long, if differences between the various experimental and numerical results are to be reconciled.

The long-time simulation of the flow revealed that, past the stage of linear instability growth, low-frequency unsteadiness of the triple point and consequently of the separation bubble sets in, having a time period of oscillation of 46 T. The long-time spatio-temporal evolution of the flow at the triple point $T_2$ of $\sim 190$T revealed the presence of spanwise corrugation of the shock as well as sinusoidal oscillations in time. The oscillation frequency corresponds to Strouhal number of $St \sim 0.028$, consistent with the existing literature on turbulent SBLI.