1 Introduction
1.1 General overview
Shockwave/turbulentboundarylayer interactions (STBLI) occur in many highspeed applications such as air intakes, control surfaces and overexpanded nozzles. If the interaction is strong enough, the flow separates, resulting in a highly unsteady flow. Of particular interest are lowfrequency components, i.e. phenomena with pertinent time scales much larger than those associated with the incoming turbulence. The associated dynamics, including unsteady thermomechanical loading and density fluctuations, can excite adverse structural responses, lead to loss of control authority, result in inefficient or even catastrophic inlet dynamics and exacerbate aerooptical distortion. Reviews by Dolling (Reference Dolling2001) and more recently Gaitonde (Reference Gaitonde2015) discuss progress in STBLI research and the challenges that lie ahead. Both emphasize lowfrequency unsteady aspects as areas of major concern, requiring control techniques for their moderation. However, characterizing the dynamics of these phenomena, especially distinguishing causation from correlation, remains a challenging task – this has inhibited the development of an optimal strategy for flow control.
Recent representative experimental (Dupont, Haddad & Debiève Reference Dupont, Haddad and Debiève2006; Dupont et al. Reference Dupont, Piponniau, Sidorenko and Debiève2008; Ganapathisubramani, Clemens & Dolling Reference Ganapathisubramani, Clemens and Dolling2009; Piponniau et al. Reference Piponniau, Dussauge, Debieve and Dupont2009; Souverein et al. Reference Souverein, Dupont, Debiève, Van Oudheusden and Scarano2010) and computational (Pirozzoli & Grasso Reference Pirozzoli and Grasso2006; Touber & Sandham Reference Touber and Sandham2009; Agostini et al. Reference Agostini, Larchevêque, Dupont, Debiève and Dussauge2012; Grilli et al. Reference Grilli, Schmid, Hickel and Adams2012; Priebe & Martín Reference Priebe and Martín2012; Aubard, Gloerfelt & Robinet Reference Aubard, Gloerfelt and Robinet2013; Morgan et al. Reference Morgan, Duraisamy, Nguyen, Kawai and Lele2013; Priebe et al. Reference Priebe, Tu, Rowley and Martín2016) efforts have explored various aspects of STBLIs and several physical mechanisms have been proposed to explain the observed unsteadiness. These may be broadly classified into two main categories. The first underscores the upstream influence of the incoming turbulent boundary layer. For instance, Ganapathisubramani et al. (Reference Ganapathisubramani, Clemens and Dolling2009) demonstrated a correlation between the lowfrequency response of the surrogate separation point and largescale, low and highspeed regions in the incoming turbulent boundary layer. The second considers downstream influence via coupling between the boundary layer, recirculation region and shock. In particular, Piponniau et al. (Reference Piponniau, Dussauge, Debieve and Dupont2009) built upon the analysis of Dupont et al. (Reference Dupont, Haddad and Debiève2006) to propose a model which explains breathing of the separation bubble and lowfrequency shock motion in terms of fluid entrainment in the mixing layer, whereby fluid from the separation bubble is continuously entrained in the mixing layer, shed downstream and must be replenished at a time scale corresponding to the lowfrequency shock oscillations. This second class of mechanisms has been demonstrated to exist even when the boundary layer is statistically free from noise near the low frequency of interest (Touber & Sandham Reference Touber and Sandham2009). It is generally accepted that while both classes of mechanisms contribute to the lowfrequency shock dynamics, the second class becomes more dominant as the size of the timemean recirculation region increases (Souverein et al. Reference Souverein, Dupont, Debiève, Van Oudheusden and Scarano2010; Clemens & Narayanaswamy Reference Clemens and Narayanaswamy2014). Morgan et al. (Reference Morgan, Duraisamy, Nguyen, Kawai and Lele2013) have recently evaluated the validity of many of these mechanisms through examination of an extensive largeeddy simulation (LES) database. Although these and subsequent studies find a correlation between the size of the ‘turbulent separation bubble’ and the longitudinal shock location, the driving mechanisms and a causal description of the lowfrequency dynamics of the separation bubble remain unclear.
Additional theoretical and numerical tools, including stability analyses, have also been deployed to understand the physical origin of lowfrequency shock unsteadiness. Concerning laminar interactions, Robinet (Reference Robinet2007) found a threedimensional global instability through BiGlobal linearstability analysis that led to selfsustained lowfrequency oscillations when the incident shock exceeded a critical angle. Recently, Guiho, Alizard & Robinet (Reference Guiho, Alizard and Robinet2016) demonstrated that laminar SBLIs are globally stable to twodimensional perturbations over a wide parameter space, behaving as globally stable selective amplifiers as opposed to unstable oscillators. Sansica, Sandham & Hu (Reference Sansica, Sandham and Hu2014) investigated the lowfrequency response of a laminar SBLI to whitenoise forcing through the solution to a nonmodal initial value problem (IVP), rather than a modal eigenvalue problem (EVP) (as distinguished by Theofilis (Reference Theofilis2011)), exciting a midfrequency Kelvin–Helmholtz (K–H) response as well as a lowfrequency broadband response. This work was a natural extension of an earlier effort by Touber & Sandham (Reference Touber and Sandham2009), who perturbed the timemean flow obtained from an LES of a turbulent SBLI using white noise to find an unstable, twodimensional, global mode with a growth time scale similar to that observed in lowfrequency bubble breathing. IVP perturbation analyses of the timemean turbulent flow were also conducted by Pirozzoli et al. (Reference Pirozzoli, Larsson, Nichols, Bernardini, Morgan and Lele2010), who observed a similar unstable global mode. They corroborated their findings with a BiGlobal linearstability analysis of the timemean turbulent flow, which, in addition to an unstable nonoscillatory (zerofrequency) global mode, identified several weakly damped oscillatory modes resembling bubble breathing extracted from lowpass filtered LES flow fields. This analysis was extended by Nichols et al. (Reference Nichols, Larsson, Bernardini and Pirozzoli2017) who describe the STBLI as a weakly damped oscillator, sustained by forcing, with a nonoscillatory unstable global mode. Notably, Sartor et al. (Reference Sartor, Mettot, Bur and Sipp2015) did not observe this nonoscillatory unstable global mode in the analysis of a transonic STBLI, and they found the least stable eigenvalues unrelatable to observed experimental unsteadiness. These observations are reconcilable with the previous studies, since the Reynolds averaged Navier–Stokes (RANS) equations, not the Navier–Stokes equations, were linearized about the turbulent mean, itself obtained from a RANS calculation that converged to a steady solution. As Nichols et al. (Reference Nichols, Larsson, Bernardini and Pirozzoli2017) identified: the linearization of the RANS equations for such a solution should be globally stable. However, the same cannot be said for linearizing the RANS equations about a steady base flow that does not satisfy the steady RANS equations, which may give rise to unstable global modes (Sartor, Mettot & Sipp Reference Sartor, Mettot and Sipp2014). A unified theoretical basis for linearizing propagation equations about a steady base flow that is not a steady solution to the given propagation equations remains to be developed.
1.2 Current contribution
The present study seeks to build on this understanding of the steady, laminar or turbulentmean, base flow (also referred to as the basic state) by applying, for the first time, a perturbation analysis to the evolving, unsteady STBLI. The goal is to examine the dynamic linear response (DLR) of the flow to identify mechanisms that result in the observed lowfrequency dynamics. Specifically, we employ concepts developed for the analysis of nonlinear dynamical systems (Khalil Reference Khalil1996), including those of Lyapunov stability and Lyapunov exponents (Skokos Reference Skokos, Souchay and Dvorak2010; Pikovsky & Politi Reference Pikovsky and Politi2016). The analysis is accomplished through calculation of the evolution of forced and constrained linear perturbations of the timeresolved STBLI base flow. This technique for analysing flow stability may be classified as global (Theofilis Reference Theofilis2011), as opposed to local (Huerre & Monkewitz Reference Huerre and Monkewitz1990), and nonmodal (Schmid Reference Schmid2007), as opposed to modal (Schmid & Henningson Reference Schmid and Henningson2001), following terminology described in the aforementioned references; it considers the global influence of the flow, and no ansatz of frequencyisolated modes is assumed, since intermodal frequency mixing will generally occur for linear perturbations of an unsteady base flow; such an occurrence may arise without any nonlinear interactions.
1.3 Dynamic linear response
In this framework, the base flow is the evolving, unsteady STBLI, as obtained from an LES of the Navier–Stokes equations in suitable fashion; the approach employed in this work is discussed in § 2. The dynamic linear response is then the statisticalensemble description, including timelocal, timemean and spectral components, of the linearperturbations evolving around the unsteady base flow subject to forcing; considerations concerning the extraction of coherent features of the turbulent flow and ensemble averaging are addressed by Hussain (Reference Hussain1986) and Jeong et al. (Reference Jeong, Hussain, Schoppa and Kim1997), respectively. Details of all aspects of the technique employed are presented in § 3. The dynamic linear response (§ 3.1) provides insight into the tendency of the evolving flow at different points in the lowfrequency cycle, in turn aiding in the determination of causal mechanisms in the STBLI. Furthermore, through proper processing and interpretation of perturbation evolution, correlations between features of unsteadiness in the low, mid and highfrequency ranges can be discerned, as can the manner in which these linear tendencies contribute to the restoration (or otherwise) of the flow to its timemean state.
Requirements for admissibility of perturbations for the dynamic linear response are similar to those for traditional stability analyses; they must satisfy the linearized Navier–Stokes (LNS) equations subject to forcing with specified boundary and initial conditions. With regard to the form of forcing, stability literature considers time harmonic, impulsive and stochastic variants (Schmid Reference Schmid2007). Our strategies are discussed in § 3.2 and include both bandlimited white noise as well as ‘native’ forcing, which has the same spectral form as the local turbulent fluctuations. A key distinction of the present effort from traditional analyses is that since the base flow of interest is highly unsteady, the linearization must be performed about the timeevolving turbulent flow. The numerical technique used solves the LNS equations through a modification of the synchronized LES procedure (SLES) described by Unnikrishnan & Gaitonde (Reference Unnikrishnan and Gaitonde2015, Reference Unnikrishnan and Gaitonde2016), who developed it to track the evolution of wave packets originating from turbulent fluctuations in the core of a jet, to understand intermittent events in the acoustic near field. Like other IVP perturbation methods, results from SLES must be processed as a statistical ensemble and do not directly reveal modal information. However, given the prohibitive computational expense of methods which analyse threedimensional steady base flows, e.g. TriGlobal (EVP) analyses (Theofilis Reference Theofilis2011), IVP methods are perhaps the only suitable techniques by which to probe timedependent threedimensional base flows.
Stability studies of timedependent flows are relatively rare, and are generally focused on the nonmodal stability of laminar flows. Timedependent nonmodal approaches may be generally classified as either examining the fundamental solution operator (Schmid Reference Schmid2007) or the adjoint equations, which are reviewed by Luchini & Bottaro (Reference Luchini and Bottaro2014). The success of these approaches in understanding timedependent laminar base flows is discussed by Schmid (Reference Schmid2007). Briefly, the dynamic linear response facilitates exploration of the nature of the fundamental solution operator using an ensemble of IVPs, each corresponding to a specific realization of the base flow and forcing. Generally, the fundamental solution operator can describe the nonmodal evolution of perturbations of steady or unsteady base flows, and it is often referred to as the Green’s function (Luchini & Bottaro Reference Luchini and Bottaro2014) or the propagator (Farrell & Ioannou Reference Farrell and Ioannou1996); it is analogous to a timedependent resolvent operator (McKeon & Sharma Reference McKeon and Sharma2010) but is not restricted to harmonic behaviour in time. Specifically, for the complex flow considered, the fundamental solution operator is not explicitly calculated, as this would be prohibitive in an unsteady calculation with $\mathit{O}(10^{8})$ discrete degrees of freedom. Rather, the linear response is determined from a statisticalensemble analysis of the propagated linear perturbations. That is, the linear response is identified through analysis of many realizations of linear perturbations (the evolution of which is specific to each forcing and baseflow realization) by invoking flow ergodicity over long duration.
The application of synchronized LES analysis to the wallbounded, turbulent, separated flow encountered in STBLIs is not trivial since there are relatively large reversed flow regions and thin, highgradient, subsonic layers near the wall, the chaotic nature of which leads to large perturbation growth rates. Perturbations growing to nonlinear magnitudes indicate that the synchronized simulations (the perturbed and baseflow realizations of the STBLI) decorrelate in time, as discussed in § 3.3. While these realizations remain statistically similar, the instantaneous difference between the two is no longer representative of the dynamic linear response of the STBLI. Since the time scales of perturbation growth to nonlinear magnitudes are later shown to be much shorter than the time scales of lowfrequency mechanisms, a difficulty arises in acquiring the dynamic linear response of lowfrequency mechanisms associated with the separation bubble and concomitant shock motion. To resolve this issue, § 3.3 introduces a procedure to linearly constrain perturbation growth and maintain linearity with respect to the base flow, enabling sustained longtime acquisition of the desired linearperturbation data. This corresponds to the implicit addition of a forcing term to the governing equations that highlights the timelocal tendency toward rapid linear growth, i.e. the influence of faster growing perturbations persists in time, while the influence of slower growing perturbations attenuates. The overall procedure of using two LES to extract the dynamic linear response is presented in § 3.4. When the linear constraint is applied, the SLES is then designated LCSLES. In the context of Lyapunov stability, the effect of this linear constraint amounts to Lyapunov renormalization, which facilitates the identification of unstable linear dynamics associated with the largest Lyapunov exponent of the flow, as discussed in § 3.5 (Skokos Reference Skokos, Souchay and Dvorak2010; Pikovsky & Politi Reference Pikovsky and Politi2016).
1.4 Application to STBLI
The linearly constrained analysis is applied to the impinging STBLI, with emphasis on understanding upstream influence, separation bubble dynamics, vortex shedding and lowfrequency shock motion. The specific problem considered is described in § 4, which includes an assessment of the base flow. Briefly, we consider an oblique shock wave generated by a $9^{\circ }$ wedge impinging on a Mach 2.33 (momentum thickness Reynolds number: $Re_{\unicode[STIX]{x1D703}}=2300$ ) turbulent boundary layer, previously studied computationally by Mullenix & Gaitonde (Reference Mullenix and Gaitonde2013) and validated by comparison with experimental data obtained at The Ohio State University Gas Dynamics and Turbulence Laboratory (Webb, Clifford & Samimy Reference Webb, Clifford and Samimy2013). The effect of stronger interaction strength is also addressed with an $11^{\circ }$ wedge shock.
The results in § 5 first consider preliminary aspects of constrained and unconstrained perturbation evolution and subsequently address the physics of the lowfrequency oscillations. For this, several cases are delineated in § 5.1, differing from each other primarily in the nature, location and duration of forcing. In § 5.2, the transient amplification without constraint is discussed to illustrate sensitivity of the chaotic flow to initial conditions, the extent of the linear growth domain, and to motivate and justify the constraining parameters chosen for subsequent work. Section 5.3 describes the statistically stationary state of the linearly constrained perturbations. The path by which the stationary state is achieved is first discussed in § 5.3.1, followed by a demonstration that the constraining technique recovers linearity in § 5.3.2. The response to different forcing parameters reveals the presence of selfsustaining linearly constrained perturbations in § 5.3.3. These observations are connected to the concepts of absolute and convective global instabilities discussed by Huerre & Monkewitz (Reference Huerre and Monkewitz1990) and Chomaz (Reference Chomaz2005), which we observe in the STBLI and undisturbed turbulent boundary layer, respectively; in response to forcing, absolutely unstable flows exhibit selfsustained oscillations, whereas convectively unstable flows respond as selective amplifiers. Additionally, these findings are related to the meanflow stability results of Touber & Sandham (Reference Touber and Sandham2009), Pirozzoli et al. (Reference Pirozzoli, Larsson, Nichols, Bernardini, Morgan and Lele2010) and Nichols et al. (Reference Nichols, Larsson, Bernardini and Pirozzoli2017), in the context of shock motion asymmetry observed by Piponniau et al. (Reference Piponniau, Dussauge, Debieve and Dupont2009). Again employing the terminology of stability theory, the ‘receptivity’ of the turbulent flow is characterized in § 5.3.4 by examining the effect of forcing location on the perturbation field, using wallnormal velocity as a sample variable; the insensitivity to forcing facilitates confirmation of the absolute nature of the instability. We examine and correlate the spectral content of both the baseflow and perturbation fields in § 5.3.5, including the attenuation factor, which linearly enforces the perturbation constraint.
The final investigation in § 5.3.6 simultaneously decomposes the turbulent baseflow and perturbation fields through bandpass temporal filtering to generate insight into mechanisms that sustain coherent motion and lowfrequency unsteadiness. Several bands are considered, each motivated by prominent frequency ranges identified in the literature, which collapse reasonably well when described by a Strouhal number based on the nominal length of the separation bubble, $L$ : $St_{L}=fL/U_{\infty }$ (Dussauge, Dupont & Debiève Reference Dussauge, Dupont and Debiève2006). Commonly, the lowfrequency dynamics $(St_{L}\sim 0.03)$ are of most interest as this band corresponds to the relatively largestscale coherent motion of the reflected shock. High–midfrequency Kelvin–Helmholtz shedding along the separated shear layer $(St_{L}\sim 0.5)$ has also been observed and discussed (Agostini et al. Reference Agostini, Larchevêque, Dupont, Debiève and Dussauge2012; Aubard et al. Reference Aubard, Gloerfelt and Robinet2013). Additionally, we identify important dynamics at a low–mid frequency $(St_{L}\sim 0.1)$ , corresponding to oscillations of the separation bubble. In fact, $St_{L}\sim 0.1$ content is found in the recirculation region power spectra of many studies (Dupont et al. Reference Dupont, Haddad and Debiève2006; Touber & Sandham Reference Touber and Sandham2009; Aubard et al. Reference Aubard, Gloerfelt and Robinet2013), but it is often not discussed in detail. The highfrequency $(St_{L}\gtrsim 1)$ response of the shock to finescale boundary layer turbulence results in the least coherent (jittering) motion and is not considered a primary concern for analysis or control efforts.
We also discuss in § 5.3.6 the bandisolated dynamics of the base flow and perturbations in the context of a massdepletion mechanism, e.g. the conceptual model of Piponniau et al. (Reference Piponniau, Dussauge, Debieve and Dupont2009). Such a model is consistent with the lowfrequency dynamics; however, like Morgan et al. (Reference Morgan, Duraisamy, Nguyen, Kawai and Lele2013) we find that the ‘turbulent separation bubble’ responds more significantly at higher frequencies $(St_{L}\sim 0.1)$ than predicted by the model $(St_{L}\sim 0.03)$ . We discuss these observations in conjunction with the results of Kiya & Sasaki (Reference Kiya and Sasaki1985), who observe higherfrequency $(St_{L}\sim 0.1)$ content in lowspeed separation bubbles.
Given the difficulty of trying to verify the causality of any single proposed mechanism driving the lowfrequency shock unsteadiness, several authors have applied simpler stochastic models, describing the motion in terms of a firstorder stochastic ordinary differential equation (ODE). Plotkin (Reference Plotkin1975) first postulated that the lowfrequency motion could be described, assuming that the shock displaced from its timemean position would be subject to a linear restoring mechanism. This model was later verified experimentally by Poggie & Smits (Reference Poggie and Smits2001) and more recently derived analytically with supporting assumptions validated through LES by Touber & Sandham (Reference Touber and Sandham2011). Sansica, Sandham & Hu (Reference Sansica, Sandham and Hu2016) also reproduced many of the key features of the oscillation by suitably forcing a laminar SBLI, showing that broadband lowfrequency upstream content is not necessary to elicit a lowfrequency response at the separation point. Rather the response can result from ‘forcing’ via transition downstream, which is then lowpass filtered by the interaction. To assess the accuracy of the assumptions employed in development of simplified models for lowfrequency unsteadiness, § 5.3.6 also examines the perturbation dynamics at different phases of the lowfrequency shock motion cycle to identify linear tendencies for restoration toward or departure from the mean state. We discuss the accuracy of the assumptions of the ODE model with these bandisolated linear tendencies, identifying phases of the lowfrequency motion that exhibit behaviour consistent or inconsistent with the model. The paper concludes by summarizing the key results in § 6.
2 Theoretical and numerical model
The compressible Navier–Stokes equations are solved in strong curvilinear nondimensional form using an LES method discussed extensively by Garmann (Reference Garmann2013). The scales used for nondimensionalization include the free stream density, $\tilde{\unicode[STIX]{x1D70C}}_{\infty }$ , the free stream velocity, $\tilde{U} _{\infty }$ , the reference boundary layer thickness, $\tilde{\unicode[STIX]{x1D6FF}}_{0}$ , the free stream temperature, $\tilde{T}_{\infty }$ , and the free stream molecular viscosity, $\tilde{\unicode[STIX]{x1D707}}_{\infty }$ , which is obtained from Sutherland’s law, $\unicode[STIX]{x1D707}=T^{3/2}(1+C_{1}/T+C_{1})$ , where $C_{1}=110.56~\text{K}/\tilde{T}_{\infty }$ is the nondimensional Sutherland’s constant and dimensional quantities are denoted by a tilde. The nondimensional variables (time, velocity, density, pressure, temperature and viscosity) are then:
The main nondimensional parameters are the Reynolds number, $Re_{\infty }=\tilde{\unicode[STIX]{x1D70C}}_{\infty }\tilde{U} _{\infty }\tilde{\unicode[STIX]{x1D6FF}}_{0}/\tilde{\unicode[STIX]{x1D707}}_{\infty }$ , and the Mach number, $M_{\infty }=\tilde{U} _{\infty }/\sqrt{\unicode[STIX]{x1D6FE}\tilde{p}_{\infty }/\tilde{\unicode[STIX]{x1D70C}}_{\infty }}$ . The ratio of specific heats is taken as constant, $\unicode[STIX]{x1D6FE}=1.4$ , and the Prandtl number, $Pr=\tilde{\unicode[STIX]{x1D707}}\tilde{C}_{p}/\tilde{k}$ , is assumed to be 0.72, where $\tilde{C}_{p}$ is the specific heat capacity at constant pressure and $\tilde{k}$ is the thermal conductivity. The thermodynamic variables are related through the perfect gas equation, $p=\unicode[STIX]{x1D70C}T/\unicode[STIX]{x1D6FE}M_{\infty }^{2}$ .
The governing equations are transformed into curvilinear coordinates, $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D709}(x,y,z)$ , $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}(x,y,z)$ and $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}(x,y,z)$ , with a transformation Jacobian, $\text{J}=\unicode[STIX]{x2202}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702},\unicode[STIX]{x1D701},\unicode[STIX]{x1D70F})/$ $\unicode[STIX]{x2202}(x,y,z,t)$ , and can be expressed in strong curvilinear nondimensional form
where $\unicode[STIX]{x1D731}$ is the conserved variable vector, $\unicode[STIX]{x1D731}=[\unicode[STIX]{x1D70C},\unicode[STIX]{x1D70C}u,\unicode[STIX]{x1D70C}v,\unicode[STIX]{x1D70C}w,\unicode[STIX]{x1D70C}E]^{\text{T}}$ , $\boldsymbol{f}$ etc. are the inviscid flux vectors, $\boldsymbol{f}^{\boldsymbol{v}}$ etc. are the viscous flux vectors and $\boldsymbol{Q}$ is the source vector. The flux vectors can be expressed as
with contravariant velocity component, $U=\unicode[STIX]{x1D709}_{t}+\unicode[STIX]{x1D709}_{x}u+\unicode[STIX]{x1D709}_{y}v+\unicode[STIX]{x1D709}_{z}w$ , and specific energy density
The deviatoric stress tensor and heat flux vector are given by
and
where Stokes’ hypothesis is assumed for the bulk viscosity coefficient, $\unicode[STIX]{x1D706}=2/3\unicode[STIX]{x1D707}$ .
Time is discretized using a variant of the implicit, approximately factored, secondorder method of Beam & Warming (Reference Beam and Warming1978) in the diagonalized form of Pulliam & Chaussee (Reference Pulliam and Chaussee1981). In this work a nondimensional time step of $\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}=0.001$ is used, and three Newtonlike subiterations are performed in the implicit scheme, at each time step, to recover accuracy due to errors introduced by linearization, factorization and explicit updating of boundary conditions. The inviscid Roe fluxes are calculated from a fifthorder, weighted essentially nonoscillatory (WENO) reconstruction of characteristic variables. This WENO scheme includes a downwind candidate stencil in the reconstruction to reduce the amount of numerical dissipation introduced in shockfree regions (Pirozzoli Reference Pirozzoli2011), compared to the original, fully upwinded, WENO scheme (Jiang & Shu Reference Jiang and Shu1996), and improve the bandwidth resolution relative to the original scheme, while maintaining fifthorder formal accuracy. The scheme degenerates to a monotonic upwind scheme for conservation laws (MUSCL) reconstruction at the second and first points from the block boundaries; however, at block interfaces, a sufficient overlap is included to maintain high order. Further details of the WENO scheme are discussed by Mullenix & Gaitonde (Reference Mullenix and Gaitonde2011). The viscous fluxes are calculated using a sixthorder, compactcentraldifference, spectrallike scheme (Lele Reference Lele1992; Visbal & Gaitonde Reference Visbal and Gaitonde2002).
No explicit subgridscale (SGS) turbulence model is employed. The rationale for this is based on several points. The efforts of Kawai, Shankar & Lele (Reference Kawai, Shankar and Lele2010), show that for wellresolved flows at low to moderate Reynolds number, such as the flow under consideration, the addition of an SGS model is not beneficial, as it introduces excessive artificial dissipation to the resolved turbulence. Similar observations for lowspeed stalled airfoils are presented by Garmann, Visbal & Orkwis (Reference Garmann, Visbal and Orkwis2013). Indeed, Priebe & Martín (Reference Priebe and Martín2012) designate results with a similar highorder WENObased approach as a direct numerical simulation (DNS). As noted by Spalart (Reference Spalart2000), when using implicit LES approaches, it is necessary to ensure that the nearwall region is sufficiently resolved, a condition that is facilitated by the choice of the low Reynolds number and highorder numerical scheme.
3 Constrained linearization technique
Linearly constrained synchronized largeeddy simulation (LCSLES) analysis can most appropriately be interpreted as a constrained linearization about the unsteady turbulent base flow. In this section we first discuss the concept of linearization about a timedependent base flow and the modifications to the linearized equations due to forcing and constraining terms, before describing how this procedure is implemented using two synchronized LES.
3.1 Governing equations of dynamic linear response
The discretized Navier–Stokes equations may be represented as
where $\unicode[STIX]{x1D6F7}_{i}$ denotes the conserved variable corresponding to the $i$ th degree of freedom of the discretized flow, $\unicode[STIX]{x1D731}$ the complete flow state and $F_{i}$ the Navier–Stokes spatial operator applied to the $i$ th degree of freedom. The total number of degrees of freedom (time excluded) corresponds to the product of the spatial discretization count with the five conserved flow variables. We seek a linearized solution, $\unicode[STIX]{x1D731}^{\boldsymbol{L}}$ , where linearization is performed in time about the evolving LES baseflow solution, $\unicode[STIX]{x1D731}^{\boldsymbol{B}}$ , resulting in a linear perturbation from the evolving base flow, $\unicode[STIX]{x1D731}^{\boldsymbol{P}_{\boldsymbol{L}}}\equiv \unicode[STIX]{x1D731}^{\boldsymbol{L}}\unicode[STIX]{x1D731}^{\boldsymbol{B}}$ . These solution vectors can be obtained from
where the time derivative of the linearized conservative variable vector is represented by a firstorder expansion about the evolving baseflow Navier–Stokes operator, and the baseflow linearization Jacobian,
of the linearized system evolves in time nonlinearly as a function of the baseflow solution. We emphasize that this method is a generalization of IVP perturbation methods applied to the timemean flow (Touber & Sandham Reference Touber and Sandham2009; Sansica et al. Reference Sansica, Sandham and Hu2014; Waindim, Bhaumik & Gaitonde Reference Waindim, Bhaumik and Gaitonde2016) in which the baseflow linearization Jacobian is a timedependent function of the evolving base flow instead of a timeindependent function of the steady timemean flow.
As noted in the introduction, this method conceptually extends a common nonmodal technique for analysing the stability of timedependent laminar flows to the fully turbulent regime. The solution to a set of equations described by a linear partial differential operator, such as those which govern the evolution of discretized linear perturbations (3.4), can be written in terms of a fundamental solution operator. For example, considering a case with forcing, $\unicode[STIX]{x1D731}^{\boldsymbol{P}_{\boldsymbol{L}}}(t)$ can be written in terms of the fundamental solution operator, $G_{ij}(t,\unicode[STIX]{x1D70F})$ , which satisfies $[\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t{\mathcal{F}}(t)]_{ij}G_{jk}(t,\unicode[STIX]{x1D70F})\equiv \unicode[STIX]{x1D6FF}_{ik}(t\unicode[STIX]{x1D70F})$ , where $\unicode[STIX]{x1D6FF}$ is the product of Dirac and Kronecker deltas, such that
The integration can proceed with a supplied initial condition, $\unicode[STIX]{x1D731}^{\boldsymbol{P}_{\boldsymbol{L}}}(t_{0})$ , to give the complete solution if the fundamental solution operator is explicitly known in the time interval $[t_{0},t]$ :
Additionally, the fundamental solution operator satisfies $\text{d}G_{ij}(t,\unicode[STIX]{x1D70F})/\text{d}t={\mathcal{F}}_{ik}(t)G_{kj}(t,\unicode[STIX]{x1D70F})$ in the interval $[\unicode[STIX]{x1D70F},t]$ with ${\mathcal{G}}_{ij}(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D70F})=\unicode[STIX]{x1D6FF}_{ij}$ , where $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta. While the fundamental solution operator can be calculated directly and subjected to further analysis in the study of laminar flows (Schmid Reference Schmid2007), memory requirements render this impractical for timedependent turbulent flows with $\mathit{O}(10^{8})$ discrete degrees of freedom. However, through ensemble observation of $\unicode[STIX]{x1D731}^{\boldsymbol{P}_{\boldsymbol{L}}}(t)$ with a variety of unique forcing and initial conditions, we can estimate the timelocal properties of the fundamental solution operator.
We now discuss two additional necessary components of the method, deferring details to subsequent sections. First, to introduce perturbations, we adopt a forcing approach which is continuous in time, represented by $\boldsymbol{S}$ . Alternatively, one could examine a set of realizations, each with a perturbed initial condition; however, the former approach is preferred in this work as it reduces the computational complexity of analysing an unsteady flow. Second, for some cases, it is necessary to constrain the growth of perturbations to enforce $\Vert \unicode[STIX]{x1D731}^{\boldsymbol{P}_{\boldsymbol{L}}}\Vert \ll \Vert \unicode[STIX]{x1D731}^{\boldsymbol{B}}\overline{\unicode[STIX]{x1D731}^{\boldsymbol{B}}}\Vert$ for linear evolution; i.e. the perturbation magnitude remains much less than the magnitude of the baseflow fluctuations. This is accomplished through inclusion of a linear attenuation factor, $\unicode[STIX]{x1D6FC}$ . The discretized forced and linearly constrained (attenuated) perturbation field then evolves in the following way:
where $\unicode[STIX]{x1D6FF}_{ij}$ represents the Kronecker delta.
3.2 Forcing term
Two independent forms of the forcing term $\boldsymbol{S}$ are considered:
The first is a scaleddown turbulent fluctuation native to the base flow, where $\overline{\unicode[STIX]{x1D6F7}_{i}^{B}}$ is the timemean base flow, while the second is whitenoise forcing, where $\unicode[STIX]{x1D6F7}_{i\,random}$ is a random real variable distributed uniformly over $[1,1]$ . In both cases $W_{i}[\boldsymbol{x},t]\in \{0,1\}$ acts to window the source region in space and time, while $\unicode[STIX]{x1D716}_{i}[\boldsymbol{x},t]$ scales the source magnitude. Generally, within the source window, the scaling function is bounded by the target perturbation magnitude: $\unicode[STIX]{x1D716}_{i}[\boldsymbol{x},t]\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}\lesssim \unicode[STIX]{x1D6F7}_{i}^{P_{L\,target}}$ . We note that scaled native fluctuation forcing does not explicitly specify a source term for the perturbations. Rather, it specifies the effect of an implicit source term in the chosen window. Thus, where $W_{i}[\boldsymbol{x},t]\neq 0$ , $\unicode[STIX]{x1D6F7}_{i}^{P_{L}}=\unicode[STIX]{x1D716}_{i}[\boldsymbol{x},t](\unicode[STIX]{x1D6F7}_{i}^{B}\overline{\unicode[STIX]{x1D6F7}_{i}^{B}})$ ; i.e. the turbulent fluctuation native to the base flow is specified as the perturbation. The results of Unnikrishnan & Gaitonde (Reference Unnikrishnan and Gaitonde2016) indicate that transients associated with the rampup of forcing (whether infinitesimal or finite duration) wash out relatively quickly. Likewise, the results are not sensitive to the details of the edges of the forcing regions, which can be sharp or diffuse.
3.3 Constraining term
Perturbations can be linearly constrained by including a real attenuation factor, $\unicode[STIX]{x1D6FC}\in \mathbb{R}$ , which operates linearly on all system degrees of freedom to damp perturbation growth. Mathematically, this reduction of perturbation growth rate coincides with lowering the real component of the eigenspectrum of the baseflow linearization Jacobian (3.5) by $\unicode[STIX]{x1D6FC}$ . Attenuation allows for the statistical analysis of linear perturbations in the unsteady environment which would otherwise grow too rapidly to be observed over a statistically significant duration of time. Attenuation only permits the effects of the most rapidly growing perturbations to persist in the unsteady flow; all weakly growing perturbations are damped out. Several strategies for attenuation have been examined and are distilled below.
The linear attenuation factor, $\unicode[STIX]{x1D6FC}$ , can take one of three general forms:
The rationale for each is as follows.

(i) No attenuation of the perturbation field in time: the perturbations evolve linearly, and their absolute growth rate is unaltered. This is the most natural case and has been the method of choice for all simulations to date, which have focused on free shear layers. The approach has the advantages of simplicity and ease of interpretation, but post facto confirmation must be exercised to ensure that $\unicode[STIX]{x1D716}_{i}$ is small enough to yield $\Vert \unicode[STIX]{x1D731}^{\boldsymbol{P}_{\boldsymbol{L}}}\Vert \ll \Vert \unicode[STIX]{x1D731}^{\boldsymbol{B}}\overline{\unicode[STIX]{x1D731}^{B}}\Vert$ so it is reasonable to interpret the perturbations as linear with respect to the base flow. The success of this approach for supersonic jets is documented by Unnikrishnan & Gaitonde (Reference Unnikrishnan and Gaitonde2015, Reference Unnikrishnan and Gaitonde2016). For wallbounded turbulence, however, this method is unsuitable in obtaining a statistically stationary linearperturbation field. Due to high sensitivity of the turbulent boundary layer to perturbations, the perturbed flow and base flow eventually decorrelate as they evolve in time with perturbations growing to nonlinear magnitudes. However, as discussed in § 5.2, cases without attenuation can still be used to examine nonstationary behaviour and shorttime perturbation growth.

(ii) Constant attenuation of the perturbation field in time: in this case, the perturbations evolve linearly, their relative growth rate is unaltered, but their absolute growth rate is changed through continuous constant attenuation; i.e. the absolute growth rate is reduced through attenuation of perturbations, but the same linear attenuation is applied to each degree of freedom, so the relative growth rate of perturbations, compared with other degrees of freedom, is preserved. Attenuation is applied linearly to the conserved quantities across all degrees of freedom, so the procedure effectively represents volume forcing to the linearized Navier–Stokes equations. This method is advantageous as no nonlinearities are introduced through attenuation. However, this approach can also suffer from decorrelation between the perturbed flow and base flow, albeit less dramatically than with no attenuation. Furthermore, the alteration of absolute growth rate requires special care in interpretation.

(iii) Variable attenuation of the perturbation field in time: this approach is introduced in the current work to allow the perturbations to evolve pseudolinearly in time, with $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}(t)=\unicode[STIX]{x1D6FC}(\Vert \unicode[STIX]{x1D731}^{\boldsymbol{P}_{\boldsymbol{L}}}\Vert )$ , i.e. $\unicode[STIX]{x1D6FC}$ varies in time, dependent on the perturbation field magnitude, introducing a weak nonlinearity; this nonlinearity is negligible as discussed in § 5.3.2. The relative growth rate of perturbations is unaltered, but their absolute growth rate is modulated. This method is advantageous as the perturbation field can be evolved to a statistically stationary state and statistics gathered over a long period of time without the potential for decorrelation since the linearity constraint, $\Vert \unicode[STIX]{x1D731}^{\boldsymbol{P}_{\boldsymbol{L}}}\Vert \ll \Vert \unicode[STIX]{x1D731}^{\boldsymbol{B}}\overline{\unicode[STIX]{x1D731}^{\boldsymbol{B}}}\Vert$ , is continuously enforced. Statistically stationary perturbation fields in this work are obtained with the variable attenuation factor taken as a simple function of the global minimum ratio (across all degrees of freedom) of the specific target perturbation magnitudes $(\unicode[STIX]{x1D6F7}_{i}^{P_{L\,target}})$ and the linear perturbations $(\unicode[STIX]{x1D6F7}_{i}^{P_{L}})$ as described in (3.10), where $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P_{L}}$ indicates the maximum norm of the perturbation field across all degrees of freedom and $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P_{L\,target}}$ is the (general) perturbation target magnitude. This attenuation method fixes the maximum norm of the linearperturbation field $(\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P})$ near to the target perturbation magnitude $(\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P_{L\,target}})$ . For the method to be pseudolinear, normalized results should be independent of the magnitude of $\unicode[STIX]{x1D6F7}_{i}^{P_{L\,target}}$ , provided $\unicode[STIX]{x1D6F7}_{i}^{P_{L\,target}}$ is small. Experience shows that $\unicode[STIX]{x1D6F7}_{i}^{P_{L\,target}}\in [1\times 10^{7},1\times 10^{5}]$ , nondimensionalized as in § 2, is an optimal range for ensuring pseudolinearity while maximizing decimal precision; i.e. the perturbations are constrained, such that they remain five to seven orders of magnitude smaller than the free stream conservative quantities. As with constant attenuation, the alteration of absolute growth rate requires special care in interpretation.
3.4 Synchronized LES
To solve (3.8) as posed would require recomputation of the baseflow linearization Jacobian at each time step and the development of a new solution technique. Furthermore, the effects of subgrid dissipation and discretization error would not generally be consistent between the linearized and baseflow solutions. These errors can rapidly compound in the chaotic base flow and can significantly limit accuracy when computing solution sensitivity (the baseflow linearization Jacobian) for highfidelity simulations. Similar concerns are discussed in detail by Vishnampet, Bodony & Freund (Reference Vishnampet, Bodony and Freund2015) in the context of adjointbased optimization. Instead, it is more practical to calculate a second twinflow LES, $\unicode[STIX]{x1D731}^{\boldsymbol{T}}$ , using the same solution technique as the baseflow LES to guarantee consistency, with the addition of forcing and constraining terms. The twinflow solution may be represented by the expansion:
When linearity is enforced as discussed above, the higherorder terms are negligible with $\unicode[STIX]{x1D731}^{\boldsymbol{T}}\approx \unicode[STIX]{x1D731}^{\boldsymbol{L}}$ , and the development of the method follows as above. The perturbation field, $\unicode[STIX]{x1D731}^{\boldsymbol{P}}\equiv \unicode[STIX]{x1D731}^{\boldsymbol{T}}\unicode[STIX]{x1D731}^{\boldsymbol{B}}$ , is then calculated as the difference between the twinflow and baseflow solutions as they are advanced synchronously in time. The perturbations are thus effectively propagated using the same highfidelity numerical scheme employed for the baseflow LES.
The synchronized LES solutions evolve in the following way:
The method ultimately requires (3.2) to be solved synchronously with (3.13), yielding the evolution of the base flow $(\unicode[STIX]{x1D731}^{\boldsymbol{B}})$ and perturbations $(\unicode[STIX]{x1D731}^{\boldsymbol{P}})$ in time. With the computation performed using nondimensional variables, we have $\Vert \unicode[STIX]{x1D731}^{\boldsymbol{B}}\Vert \sim \mathit{O}(1)$ . Applying variable attenuation with $\unicode[STIX]{x1D6F7}_{i}^{P_{L\,target}}=10^{6}$ , we have $\Vert \unicode[STIX]{x1D731}^{\boldsymbol{P}}\Vert \sim \mathit{O}(10^{6})$ and $\Vert \unicode[STIX]{x1D731}^{\boldsymbol{P}}\Vert ^{2}\sim \mathit{O}(10^{12})$ with $\Vert \unicode[STIX]{x1D731}^{\boldsymbol{B}}\Vert \gg \Vert \unicode[STIX]{x1D731}^{\boldsymbol{P}}\Vert \gg \Vert \unicode[STIX]{x1D731}^{\boldsymbol{P}}\Vert ^{2}$ , justifying a linearized interpretation over approximately six digits of precision. This linearly constrained synchronized largeeddy simulation (LCSLES) method with variable attenuation provides a reasonable approximation to the solution linearized and constrained around the evolving base flow, in which discretization error is consistent between the linearized and baseflow solutions. This allows for the study of perturbation evolution in fully turbulent environments posed in terms of a nonmodal initial value problem.
3.5 Relationship between constrained and unconstrained perturbations
To aid in understanding and interpreting constrained perturbations, we briefly discuss their relation to unconstrained perturbations in the context of Lyapunov stability (Khalil Reference Khalil1996; Skokos Reference Skokos, Souchay and Dvorak2010; Pikovsky & Politi Reference Pikovsky and Politi2016). Consider the evolution, without forcing, of an initial unconstrained perturbation impulse, $\unicode[STIX]{x1D731}^{\boldsymbol{u}}(t_{0})$ , satisfying $\Vert \unicode[STIX]{x1D731}^{\boldsymbol{u}}(t_{0})\Vert _{\infty }=\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P_{L\,target}}$ . We may define a Lyapunov norm for the fundamental solution operator that propagates this impulse, unconstrained, in time, as the ratio of maxnorms of the perturbation at a future time relative to the initial perturbation impulse:
The last equality includes the integrated attenuation factor, $\unicode[STIX]{x1D6FC}(t)$ , corresponding to the same baseflow realization and initial perturbation impulse, but evolving subject to constraint, with target perturbation magnitude $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P_{L\,target}}$ . This equality holds, noting that $\unicode[STIX]{x1D6FC}(t)$ is defined to enforce $\Vert \unicode[STIX]{x1D731}^{\boldsymbol{c}}(t)\Vert _{\infty }/\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P_{L\,target}}=1$ ; i.e. the gain of the constrained system is unity, and the gain of the corresponding unconstrained system (the Lyapunov norm) grows exponentially with the timeintegrated attenuation factor. Incidentally, the largest Lyapunov exponent of the unconstrained system corresponding to this initial impulse is
The fundamental solution operator of the constrained linear system can then be related to that of the corresponding unconstrained system,
so that the propagation of constrained perturbations can be described:
where $\unicode[STIX]{x1D731}^{\boldsymbol{u}}(t,t_{0})$ describes the initial impulse of the unconstrained system that would result in the same perturbation state $(\unicode[STIX]{x1D731}^{\boldsymbol{u}}(t)=\unicode[STIX]{x1D731}^{\boldsymbol{c}}(t))$ at time, $t$ , as the initial impulse of $\unicode[STIX]{x1D731}^{\boldsymbol{c}}(t_{0})$ in the constrained system. Notably, $\Vert \unicode[STIX]{x1D731}^{\boldsymbol{u}}(t,t_{0})\Vert _{\infty }=\exp \int _{t_{0}}^{t}\!\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D70F})\,\text{d}\unicode[STIX]{x1D70F}$ . Thus, the constrained linear system faithfully reproduces the dynamics of the corresponding unconstrained linear system, subject to the scaling factor of the unconstrained Lyapunov norm. The process described in this section is commonly known as Lyapunov renormalization (Wolf et al. Reference Wolf, Swift, Swinney and Vastano1985), through which the most rapidly growing perturbations of the unconstrained system (those associated with the largest Lyapunov exponent) manifest as constant magnitude perturbations in the constrained system, and weakly growing perturbations in the unconstrained system are attenuated in the constrained system. This process is a special case of repetitive Gram–Schmidt orthonormalization, through which the entire Lyapunov spectrum may be obtained (Benettin et al. Reference Benettin, Galgani, Giorgilli and Strelcyn1980; Wolf et al. Reference Wolf, Swift, Swinney and Vastano1985; Geist, Parlitz & Lauterborn Reference Geist, Parlitz and Lauterborn1990; Skokos Reference Skokos, Souchay and Dvorak2010).
The dynamic linear response (DLR), which elicits information about the linear stability of a turbulent flow, or any timedependent nonlinear dynamical system, may draw connotative association with the common technique termed dynamic mode decomposition (DMD) (Rowley et al. Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009; Schmid Reference Schmid2010); however, the two are quite different methods. The former concerns the propagation of linear perturbations about the unsteady turbulent flow and is necessarily an in situ analysis, while the latter is a databased postprocessing technique, which has connections to discrete Fourier modes, Koopman modes and linear global modes, in certain limiting cases. Commonly, when DMD is performed on fluctuations of a statistically stationary turbulent flow, Fourier modes result (Chen, Tu & Rowley Reference Chen, Tu and Rowley2012); when DMD is performed on a linear system with a timeindependent propagator, global modes result, and the timeindependent linear Koopman operator, not a linear projection, is approximated (Schmid Reference Schmid2010). Notably, DMD could be applied to fluctuations of the base flow to identify Fourier modes; however, as the baseflow propagator is nonlinear and timedependent, DMD modes have no realized connection to linear stability of the turbulent base flow. Therefore, the information provided in the DLR of the base flow, which can identify tendencies in the linear stability of flow trajectories, is unique from that provided by the DMD of the base flow. The application of DMD on the propagated linear perturbations will likely provide additional insight into the DLR of the flow; however, this extension for analysing the perturbations is beyond the scope of this work.
4 Description of STBLI baseflow simulation
The STBLI examined consists of an oblique shock, generated by a wedge of angle $\unicode[STIX]{x1D703}=9^{\circ }$ in the far field, impinging on a flat plate turbulent boundary layer. The flow conditions of the experiment (Webb et al. Reference Webb, Clifford and Samimy2013) and LES are summarized in table 1. For computational accessibility, the simulation Reynolds number is lowered by a factor of 10 from the experimental Reynolds number, by reducing the free stream pressure, while maintaining the free stream velocity, Mach number and boundary layer thickness. This alteration does not contaminate the analysis of desired unsteady features, as also observed in prior studies (Priebe & Martín Reference Priebe and Martín2012; Aubard et al. Reference Aubard, Gloerfelt and Robinet2013; Morgan et al. Reference Morgan, Duraisamy, Nguyen, Kawai and Lele2013).
An equilibrium turbulent boundary layer is generated through bypass transition of an incoming laminar boundary layer, which is tripped using a counterflow body force, as discussed by Waindim & Gaitonde (Reference Waindim and Gaitonde2016). The effect is introduced in the simulation as a momentum source term in (2.2). The body force is taken to be uniform in time, so as to avoid lowfrequency contamination that can introduce long lasting spurious signatures at the relatively low Reynolds numbers considered. Similarly, the force is assumed homogeneous in the spanwise direction to avoid biasing the development of any particular spanwise wavenumbers.
Concerning boundary conditions, in all cases, a noslip, zero normalpressuregradient wall is used for the wall boundary. The wall temperature is fixed at $1.95T_{\infty }$ in the region of the counterflow body force, to accelerate transition, before returning to an adiabatic state downstream, in which zero normal temperature gradient conditions are enforced, except in the simulations with spatial refinement discussed below, wherein the wall remains adiabatic throughout. The inflow is a specified Blasius boundary layer profile, which is tripped by the counterflow body force. The shock is imposed at the top boundary through direct specification of the Rankine–Hugoniot conditions. Periodicity is assumed in the spanwise direction, and all other boundary conditions are extrapolated, assuming zero boundarynormal gradients, from the interior. Grid stretching to the downstream and free stream boundaries, in combination with the nonoscillatory spatial scheme, provides for sufficient damping of fluctuations, facilitating the application of these boundary conditions without spurious reflections.
We consider several computational grids (denoted A–D) to ensure that the size of the spanwise domain as well as local grid resolution are adequate, by examining convergence of both mean as well as fluctuating quantities. Details are summarized in table 2. Noting that grid requirements depend on flow parameters as well as numerical scheme resolution, we take our cue from recent published studies. The refined grids (A and B) approach the DNS resolution of previous efforts (Priebe & Martín Reference Priebe and Martín2012), while the standard grids (C and D) are comparable to those in the STBLI LES efforts of Touber & Sandham (Reference Touber and Sandham2009), Agostini et al. (Reference Agostini, Larchevêque, Dupont, Debiève and Dussauge2012), Aubard et al. (Reference Aubard, Gloerfelt and Robinet2013). The grids employed are structured, and spacing in the interaction region is constant in the streamwise (nominal, $\unicode[STIX]{x0394}x^{+}=23.5$ ; refined, $\unicode[STIX]{x0394}x^{+}=16.4$ ) and spanwise (nominal, $\unicode[STIX]{x0394}z^{+}=8.9$ ; refined, $\unicode[STIX]{x0394}z^{+}=6.1$ ) directions. In the wallnormal direction, hyperbolic tangent growth is applied from the wall (nominal/refined, $\unicode[STIX]{x0394}y_{min}^{+}=0.4$ ) to the height of the incoming boundary layer (nominal, $\unicode[STIX]{x0394}y^{+}=4.0$ ; refined, $\unicode[STIX]{x0394}y^{+}=2.4$ ), above which a powerlaw spacing is used; the superscript ‘ $+$ ’ denotes normalization with respect to the boundary layer inner variables. To explore the effect of spanwise domain size, the width is varied between $2\unicode[STIX]{x1D6FF}_{o}$ (B and D), commonly employed in prior studies (Touber & Sandham Reference Touber and Sandham2009; Agostini et al. Reference Agostini, Larchevêque, Dupont, Debiève and Dussauge2012; Priebe & Martín Reference Priebe and Martín2012), and $5\unicode[STIX]{x1D6FF}_{0}$ (A and C). As shown below, the results indicate that a spanwise domain size of $2\unicode[STIX]{x1D6FF}_{0}$ is sufficient. The base flow and perturbations are therefore simulated primarily on grids B and D respectively. It is also evident that the nominal grid (D) adequately reproduces the unsteady characteristics. The longest duration simulations required to analyse the frequency content were therefore performed on grid D.
Characteristics of the turbulent boundary layer are presented in figure 1. Turbulence has fully developed by the streamwise region between $x/\unicode[STIX]{x1D6FF}_{0}=50$ and $x/\unicode[STIX]{x1D6FF}_{0}=65$ , which constitutes the interaction region for the STBLI. Figure 1(a) confirms agreement of the Van Driest transformed velocities with the inner law and log law at various streamwise locations. Note that a consequence of the choice of low Reynolds number is the relative absence of a pronounced wake region (Pirozzoli & Bernardini Reference Pirozzoli and Bernardini2013). This is characteristic of several other efforts in numerical simulation of STBLI (Pirozzoli & Bernardini Reference Pirozzoli and Bernardini2011; Priebe & Martín Reference Priebe and Martín2012; Aubard et al. Reference Aubard, Gloerfelt and Robinet2013; Morgan et al. Reference Morgan, Duraisamy, Nguyen, Kawai and Lele2013) and, as shown below, has relatively little influence on the features of the lowfrequency phenomena of interest. Figure 1(b) validates the normal Reynolds stresses with analytical results for the outerscale normal Reynolds stresses given by Perry & Li (Reference Perry and Li1990). Results on grid A are closer to DNS resolution (Pirozzoli, Grasso & Gatski Reference Pirozzoli, Grasso and Gatski2004), and the deviation of the streamwisenormal Reynolds stress with grids C and D is modest. Figure 1(c) presents the onedimensional energy spectra for streamwise velocity fluctuations near the end of the interaction region. A suitable equilibrium turbulent boundary layer is observed with no lowfrequency tones. The spectra further demonstrate the anticipated theoretical rolloff $\propto k^{5/3}$ at high wavenumber and plateau toward low wavenumber as anticipated for onedimensional (1D) spectra. The reader is referred to Waindim & Gaitonde (Reference Waindim and Gaitonde2016) for further details regarding properties of the incoming turbulent boundary layer.
A comparison of the timemean computational and experimental (Webb et al. Reference Webb, Clifford and Samimy2013) flow fields is presented in figure 2. Figure 2(a,b) shows timemean contours from the fine grid (B), overlaid on filled contours from the nominal grid (D), to demonstrate a reasonable degree of grid convergence, with the latter grid solution producing a slightly larger timemean separation. Corresponding particle image velocimetry (PIV) data in the same reference panel and colour map are shown in figure 2(c,d). The boundary layer thickness, interaction length ( $L_{int}$ , the streamwise distance between the wallextrapolated reflected and impinging shocks), and point of shock impingement are well matched between the PIV and LES. The experimental data have some contamination due to surface oil flow measurements and an extraneous impinging expansion fan downstream of the interaction. This expansion likely accounts for the flow acceleration and thinner postreattachment boundary layer observed in the experiment compared with the simulation; the mean flow observed in the simulation is more favourably comparable to studies without this expansion, including those by Piponniau et al. (Reference Piponniau, Dussauge, Debieve and Dupont2009) and Morgan et al. (Reference Morgan, Duraisamy, Nguyen, Kawai and Lele2013). Nonetheless, the key features of the experiment are reproduced faithfully by the simulation. The form of the separated region has a more asymmetric nature in simulation than in experiment. Similar behaviour is observed and discussed by Touber & Sandham (Reference Touber and Sandham2009), concerning a validation study that compares results at identical Reynolds numbers.
The timemean skin friction and surface pressure profiles in the interaction region are shown in figure 3(a). The interaction region is defined as that between the wallextrapolated position of the timemean reflected shock and the shock impingement point, while the separation region is that with negative timemean streamwise skin friction. The surface pressure does not show a plateau region as is suggested by Délery & Dussauge (Reference Délery and Dussauge2009), since this is a relatively weak interaction, with large unsteady shock displacements compared to separation height. Similar profiles have been reported in other LES (Touber & Sandham Reference Touber and Sandham2009) and DNS (Pirozzoli & Grasso Reference Pirozzoli and Grasso2006). The trend of skin friction is also generally consistent with the previously mentioned studies, though there is quantitative variability of reported skinfriction values in the separation region (Aubard et al. Reference Aubard, Gloerfelt and Robinet2013). Results of the current study most closely resemble the DNS of Pirozzoli & Grasso (Reference Pirozzoli and Grasso2006).
A comparison with experiment of premultiplied and normalized wall pressure power spectral density (PSD) at locations described in table 3 and figure 3(a) is shown in figure 3(b). All computational and experimental probe data exhibit lowfrequency components. The peak frequency increases moving downstream of the separation point. The most downstream probe is placed in the aft of the recirculation region, where the peak frequency is $St_{L}\sim 0.5$ , corresponding primarily to K–H dynamics. Moving upstream, the peak frequency drops and becomes more indicative of the lowfrequency $(St_{L}\sim 0.03)$ bubble breathing and shock motion. The simulations reproduce the observed trend, with the quantitative agreement becoming better for higher frequencies. Probe placement is exactly reproduced in the simulation with respect to the impinging shock location and experimental probe displacements. As such, discrepancies in comparison of point probe data are primarily due to the difference in separation length between the experiment and simulation. Overall, the results indicate clearly that the simulations capture the key unsteady phenomena observed in the experiment. Comparison with experimental schlieren in figure 3(c) also shows good agreement. Figure 3(d) illustrates the probability of reversed flow in coordinates nondimensionalized by the mean separation length ( $L_{sep}$ ), with the mean separation point ( $x_{sep}$ ) located at $x^{\ast }=0$ : $x^{\ast }=(xx_{sep})/L_{sep}$ . This metric compares well with previous simulations by Agostini et al. (Reference Agostini, Larchevêque, Dupont, Debiève and Dussauge2012) and experiments by Dupont et al. (Reference Dupont, Piponniau, Sidorenko and Debiève2008). Finally, figure 4 shows a threedimensional instantaneous rendering of the baseflow interaction. In addition to the turbulent boundary layer, which is visualized with an isosurface of the $Q$ criterion coloured by streamwise velocity, the shocks and midplane shading are indicated by an isosurface and contours of $\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p$ respectively. The midplane visualization connects the shock structure with the nearfield wave pattern which proceeds periodically downstream from the reflected shock, away from the interaction. The separated flow, the amplification of turbulence through the interaction, and largescale lowfrequency coherent displacement of the reflected shock are dominant features of the interaction.
5 Results
5.1 Forcing parameters and cases considered
Numerous simulations were performed to understand numerical and physical aspects of the dynamic linear response of the timeresolved STBLI base flow, using different forcing and constraining parameters as summarized in table 4. Forcing locations are described in terms of their absolute position in the spanwisenormal plane and are periodic in the spanwise direction, encompassing a single line of grid points. However, the forcing itself is nonuniform in the spanwise direction, varying either as a function of the base flow (native, i.e. based on local turbulent fluctuations in the base flow) or randomly (white noise). Furthermore, forcing is applied to all pertinent variables, including mass density, three components of momentum density and energy density. Thus, no bias is introduced regarding twodimensionality or obliqueness of the perturbations. As indicated in table 4, cases 1–5 rely on native forcing, while cases 6–10 use whitenoise forcing, both of which were described earlier in the context of (3.9). Previous work on noise source identification in a turbulent jet relied on the use of native forcing, since the objective was to track existing fluctuations, not to introduce foreign perturbations (Unnikrishnan & Gaitonde Reference Unnikrishnan and Gaitonde2015, Reference Unnikrishnan and Gaitonde2016). However, in this study we observe the base flow to be selective, producing very similar linearly constrained perturbation fields regardless of forcing type, so we choose also to apply random whitenoise forcing: this facilitates confirmation of the hypothesis that this selectivity is due to the base flow and not the manner of forcing.
In cases 11–20, the whitenoise forcing is turned off after a brief transient of roughly $\unicode[STIX]{x1D70F}=5$ , and the simulations are allowed to proceed; this procedure effectively introduces an impulse at $\unicode[STIX]{x1D70F}=5$ , from which the simulations proceed without forcing. These cases are designed to highlight the selfsustaining constrained perturbation field due to the presence of the STBLI: the recirculating flow admits much larger perturbation growth rates than the downstream boundary layer, due to constructive feedback through recirculation, and acts to anchor the perturbation field to the interaction region. In contrast, as noted earlier, perturbations introduced in a turbulent boundary layer, without an impinging shock, are convected out of the domain.
5.2 Analysis of transient amplification of unconstrained perturbation field
A study of flow sensitivity to different forcing magnitudes and locations with unconstrained perturbations is performed, employing native forcing with no attenuation. Several forcing locations are considered; the results are distilled with forcing in the boundary layer: $(x=53.8,y=0.79)$ (cases 1–3, which differ in forcing magnitude) and free stream: $(x=50.9,y=2.83)$ (case 4). Perturbations injected in the boundary layer tend to rapidly decorrelate the base flow and twin flow. To quantify perturbation growth, we discuss the maxnorm of the density perturbation field. Figure 5 shows that this norm grows at an exponential rate, regardless of the forcing magnitude. For forcing locations in the boundary layer, this exponential growth begins immediately and eventually plateaus after the perturbation field saturates, at which time the decorrelation between the unperturbed and perturbed states is complete; i.e. the perturbations grow to nonlinear magnitudes. For forcing locations in the free stream, the growth is relatively small until the effects of the perturbation enter the boundary layer, after which it behaves similarly to forcing in the boundary layer.
This high sensitivity to perturbation demonstrates the chaotic nature of wallbounded turbulence. Notably, the exponential nature of nonmodal perturbation growth observed for the current unsteady turbulent base flow contrasts with the common observation of nonmodal perturbation growth for steady laminar base flows. For example, Matsubara & Alfredsson (Reference Matsubara and Alfredsson2001) discuss the algebraic growth associated with streamwise elongated streaky structures in laminar boundary layers. For high levels of free stream turbulence, this nonmodal algebraic growth can occur rapidly and bypass classical modal exponential growth (Tollmien–Schlichting waves), to dominate the early stage of the transition process (Brandt, Schlatter & Henningson Reference Brandt, Schlatter and Henningson2004). However, as figure 5 demonstrates, the current nonmodal analysis of the unsteady turbulent environment suggests that exponential growth is naturally recovered in the turbulent regime, consistent with an approximately constant leading Lyapunov exponent for the dynamical system describing the timeresolved turbulent flow.
The growth rates are similar for all forcing scenarios, until the perturbation field amplitude approaches the approximate linearity threshold, above which the assumption that higherorder terms in the expansion of (3.11) are negligible becomes less valid. The linearity threshold is approximate, not definite, and depends on the base flow. For the boundary layer under consideration, nonlinear effects, which are evident in the trajectories of figure 5 as departure from exponential growth, ensue in the range $\Vert \unicode[STIX]{x1D731}^{\boldsymbol{P}}\Vert \sim \mathit{O}(10^{5},10^{3})$ ; a conservative estimate for the linearity threshold is included near $\Vert \unicode[STIX]{x1D731}^{\boldsymbol{P}}\Vert \sim \mathit{O}(10^{5})$ to approximately demarcate the linear and nonlinear ranges of perturbation growth. To ensure perturbations remain linear, the maxnorm of the perturbation field should remain much less than the magnitude of the linearity threshold. This is facilitated in § 5.3, by constraining perturbation growth, with $\unicode[STIX]{x1D6F7}_{i}^{P_{L\,target}}=10^{6}<\mathit{O}(10^{5},10^{3})$ .
Without constraint, as the perturbations grow larger than the linearity threshold, the growth rate diminishes due to nonlinearity and the perturbation field loses significance, since the twin flow and base flow are now different decorrelated realizations of the STBLI. The two realizations remain statistically similar, but the instantaneous correlation between the two tends toward zero, and the difference between the two can no longer be reasonably interpreted as a perturbation quantity. However, as with many other dynamical systems, the path to nonlinearity can yield useful insight. Information regarding local sensitivity in the boundary layer and the interaction region can be obtained by performing short duration $[\mathit{O}(\unicode[STIX]{x1D70F}\sim 1)]$ simulations which remain in the linear regime.
Finally, we note that the time scale for perturbation growth to nonlinear magnitudes in this STBLI is of order $\unicode[STIX]{x1D70F}\sim 1$ . In previous studies of a supersonic jet (Unnikrishnan & Gaitonde Reference Unnikrishnan and Gaitonde2016), nonlinearperturbation magnitudes were not encountered in simulations conducted for much longer duration, indicating the absence of strong feedback mechanisms in the turbulent jet when perturbations are injected in the supersonic flow.
5.3 Analysis of statistically stationary linearly constrained perturbation field
The relevant time scales for lowfrequency unsteadiness in the STBLI range from $\unicode[STIX]{x1D70F}\sim \mathit{O}(10)$ for Kelvin–Helmholtz phenomena, to $\unicode[STIX]{x1D70F}\sim \mathit{O}(200)$ for lowfrequency bubble breathing. Since the unconstrained perturbation field grows rapidly to a nonlinear state, to study these time scales in a statistically significant manner we use the method of § 3.3 to constrain the perturbation field to linear magnitudes. In this section we consider cases 5–20, which employ variable attenuation and ensure linearity with respect to the base flow. As discussed in § 3, the interpretation of the constrained perturbation field is not as intuitive, since it effectively results in volume forcing of the linearized Navier–Stokes equations, in a manner that highlights the most dominant growth modes.
5.3.1 Evolution of perturbation field to statistically stationary state
The evolution of the vertical velocity perturbation field (case 7) to a statistically stationary state is depicted in figure 6; the baseflow unsteady pressure field is superposed for reference. In the initial stages, the evolution is somewhat isotropic. However, once the perturbations enter the separated region, their evolution becomes more directional in the downstream direction. By panel (c), the outer envelope starts to approximate the reflected shock location. After about $10\unicode[STIX]{x1D70F}$ , the perturbations have filled the separated region, and waves with acoustic features are evident downstream. These are also strongly directional, and do not permeate upstream of the reflected shock. Figure 6(e) provides a typical representation of the instantaneous perturbation field under statistically stationary conditions. While the base flow may take many flowthrough times to reach statistical stationarity, the perturbations reach stationarity quickly, since they are a very small deviation from the base flow, and their growth is dominated by the absolute instability. The constrained perturbation field is localized to the region just upstream of the reflected shock foot and follows the characteristic along the reflected shock to the far field, demonstrating that with variable attenuation, the perturbation field protrusion upstream is limited. In contrast, for an unconstrained case, the perturbations travel everywhere, including upstream through the subsonic boundary layer, as they grow to nonlinear magnitudes. However, linearperturbation growth is much stronger in the interaction region and farther downstream than upstream, so these regions of dominant perturbation growth are highlighted when the perturbations are constrained, since weakly growing upstream perturbations are damped relative to strongly growing perturbations in the interaction region and downstream; the magnitude of the latter remain approximately constant due to variable attenuation.
Key features of the statistically stationary perturbation field include convective structures that are shed downstream from the separation region. Additional finescale perturbation structures highlight vortices present in the baseflow postreattachment boundary layer, around which sharp gradients in the flow promote perturbation growth. Also prominent are relatively large magnitude perturbations along the reflected shock wave, demonstrating the influence of perturbation dynamics of the separation region on the reflected shock. Generally, the constrained perturbation field evolves to highlight turbulent features with high potential for linear growth, away from the base flow, which is a novel way to visualize the instantaneous character of a turbulent flow.
5.3.2 Linearity of constrained perturbation field
Several of the following conclusions rely on the claim that the LCSLES procedure represents a reasonable approximation for the evolution of constrained linear perturbations as discussed throughout § 3. Notably, any method for propagating linear perturbations through a chaotic nonlinear flow incurs errors that accumulate over time, due to sensitivity of the flow to initial conditions. In this section, we quantify the magnitude of nonlinear errors introduced through the combined effects of the basic SLES procedure (i.e. calculating the evolution of linear perturbations as the difference between two evolving fully nonlinear systems), the variable attenuation procedure, and finite decimal precision: the major contributing causes of nonlinearity. For this purpose, we analyse cases 14–20 that are all initialized with a perturbation impulse taken from a statistically stationary instant of case 11. Each of cases 14–20 proceeds with a distinct perturbation target magnitude, separated by an order of magnitude as described in table 4; the initial impulse for these cases is also scaled with the target perturbation magnitude relative to case 11. Through this premise, linear perturbations would evolve identically in cases 14–20, subject to the relative scaling of the target perturbation magnitude and initial impulse; thus, nonlinear effects may be observed through divergence of the trajectories in time. These simulations evolve for a duration of $100\unicode[STIX]{x1D6FF}_{0}/U_{\infty }$ , to observe the accumulation of nonlinearity over a time scale corresponding to the order of lowfrequency unsteadiness.
The effects of nonlinearity are described quantitatively in figure 7. Figure 7(a) describes the evolution of pressure perturbations for cases 14–20, at a representative location in the outer layer of the postreattachment boundary layer; the greatest perturbation magnitudes are observed in this vicinity. The trajectories remain similar over the time window of interest, demonstrating that linearity of the method is maintained. As time progresses, finerscale features of the signal are distorted for large magnitude perturbations due to accumulation of nonlinearity, as well as for the smallest magnitude perturbation ( $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P_{L\,target}}=1\times 10^{9}$ ) due to lack of decimal precision. Figure 7(b) shows the level of accumulated nonlinear error with time for each case with respect to case 19 (this case has the smallest target perturbation magnitude for which errors due to lack of decimal precision are not evident and is therefore considered to be the most accurate linear solution). The error is lowpass filtered in time to aid in discerning trends between the level of error and the target perturbation magnitude, and the rootmeansquare (r.m.s.) level of the normalized pressure perturbations is included to provide perspective for the error magnitude. Concerning cases in which the lack of decimal precision does not introduce error, the errors initially grow from machine precision, with faster growth rates corresponding to larger perturbation target magnitudes, as anticipated. Near $\unicode[STIX]{x1D70F}\approx 35$ , an event is encountered which introduces substantial nonlinearity into all perturbations, after which the same trend in growth is recovered. Encouragingly, over the duration of these simulations, the cumulative nonlinear error for the target perturbation magnitude of $1\times 10^{6}$ (employed in cases 5, 7, 9–13, 17) remains less than ${\approx}10\,\%$ of the r.m.s. pressure perturbations. Maintaining linearity at this level of precision over such long duration in a chaotic system is not trivial, and these observations attest to the effectiveness of the LCSLES method.
5.3.3 Timemean perturbation field and global absolute linear instability
We now discuss the behaviour of the linearly constrained perturbations in the context of a global absolute linear instability. The selfsustaining nature of the perturbations (cases 11–20) indicates the existence of constructive feedback through the recirculation region. This results in large perturbation growth rates in the recirculation region that overwhelm convective growth downstream as well as any additional forcing of the interaction, as discussed further in § 5.3.4. This is consistent with the presence of a global absolute linear instability, as discussed by Chomaz (Reference Chomaz2005). Discussion in this section focuses on the timemean behaviour of this absolute instability, but it is relevant to note that the timelocal behaviour of this absolute instability is not always consistent with its timemean behaviour, as addressed in § 5.3.6.
The conclusion that selfsustaining constrained perturbations are indicative of a global absolute instability can be confirmed using the framework discussed by Huerre & Monkewitz (Reference Huerre and Monkewitz1990) and Schmid & Henningson (Reference Schmid and Henningson2001) for classifying the nature of instability based on properties of the fundamental solution operator (referred to, in these works, in a less general form: a Green’s function). Since our method is global (Theofilis Reference Theofilis2011) and nonmodal (Schmid Reference Schmid2007), including effects from the entire domain, with no restriction on perturbation form, apart from the requirement that they remain linear with respect to the expansion in (3.11), instabilities identified through this method are global and nonmodal in nature; the fundamental solution operator propagates the complete global state space (discrete) or global phase space (continuous) in time.
For timedependent base flows, we generalize the classification of instabilities as follows. Consider the discretized form of (3.7), without forcing and without constraining, which propagates linear perturbations through the flow:
A flow is classified as globally unstable if
where $\Vert \boldsymbol{{\mathcal{G}}}(t,t_{0})\Vert _{LN}$ is the Lyapunov norm maximized over all initial impulses; conceptually, this is the maximum possible gain of the system represented by (5.1). A globally unstable flow is additionally classified as absolutely unstable if
A globally unstable flow is additionally classified as convectively unstable if
where $i_{0}$ and $j_{0}$ represent timeindependent indices for the considered degree of freedom, nonboundary sets do not include $i_{0}$ on the spatial domain boundary when considering finite domains and $t_{0}$ represents the initial time of impulse perturbation to the base flow, which is arbitrary and independent of the other variables; $\Vert {\mathcal{G}}_{i_{0}j_{0}}(t,t_{0})\Vert _{\infty }$ represents the maxnorm of the fundamental solution operator as a function of the current, $t$ , and impulse, $t_{0}$ , times, over the finite set of $i_{0}$ and $j_{0}$ . Naturally, the finite, nonboundary sets of degrees of freedom discussed in (5.3)–(5.4) represent spatially localized regions of the flow. In this sense, equation (5.3) states that the flow is absolutely unstable if the location of maximum perturbation growth is confined to a single local region of the flow. Similarly, equation (5.4) states that the flow is convectively unstable if the location of maximum perturbation growth is not confined to any single local region of the flow. The above use of ‘local’ does not imply the use of any local profile in the analysis, since the full, spatially threedimensional, ‘global’ flow is analysed. In this context, when perturbation forcing is turned off after an initial impulse, selfsustaining constrained perturbations arise in the STBLI interaction region: the perturbations are perpetually growing, hence the need for constraint, which indicates that the condition, $\limsup _{t\rightarrow \infty }\Vert \boldsymbol{{\mathcal{G}}}(t,t_{0})\Vert _{LN}\rightarrow \infty$ , is satisfied for some $t_{0}$ , consistent with a global instability. Furthermore, the constrained perturbations, which reveal the most rapid nonmodal growth, remain concentrated in the interaction region, which indicates that the condition $\limsup _{t\rightarrow \infty }\Vert {\mathcal{G}}_{i_{0}j_{0}}(t,t_{0})\Vert _{\infty }/\Vert \boldsymbol{{\mathcal{G}}}(t,t_{0})\Vert _{LN}\in (0,1]$ is satisfied for some nonboundary set of $i_{0}$ , $j_{0}$ , $t_{0}$ , in which the $i_{0}$ correspond to degrees of freedom in the interaction region, consistent with a global absolute instability. (Alternatively, for a turbulent boundary layer, the constrained perturbations are not selfsustaining and convect downstream, indicating that the condition $\limsup _{t\rightarrow \infty }\Vert {\mathcal{G}}_{i_{0}j_{0}}(t,t_{0})\Vert _{\infty }/\Vert \boldsymbol{{\mathcal{G}}}(t,t_{0})\Vert _{LN}=0$ is satisfied, for all nonboundary sets of $i_{0}$ , $j_{0}$ , $t_{0}$ , consistent with a global convective instability.)
We can quantify the effect of this global absolute instability in the STBLI by examining the time mean of the selfsustaining linearly constrained perturbation field, $\overline{\unicode[STIX]{x1D731}^{\boldsymbol{P}}}$ , which we designate as the timemean linear tendency of the base flow. We interpret the nonzero timemean linear tendency as the timemean behaviour of the absolute instability in the unsteady environment. This is clarified by noting that the time mean of the perturbation field between two statistically identical turbulent realizations (twin and base) will tend toward zero with time. However, for this STBLI, the absolute instability present in the base flow manifests as a timemean linearly constrained perturbation field that does not tend toward zero with time, since the instability biases linearperturbation growth in the unsteady environment. Indeed, as noted earlier, the constraints of § 3.3 effectively highlight this absolute instability at the expense of other more weakly growing (convective) instabilities, which are damped. Since the twinflow simulation solves the forced and linearly constrained Navier–Stokes equations, it represents a linear departure from the baseflow simulation. As such, their time means (twin and base) can be slightly different – the magnitude of the difference being determined by the perturbation target magnitude ( $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P_{L\,target}}$ ) – and is reflective of the timemean behaviour of the underlying absolute instability. The underlying bias in linearperturbation growth in the unsteady environment, leading to a nonzero timemean linear tendency, can be thought of as analogous to an unstable zerofrequency (nonoscillatory) mode, which can be found in global IVP and EVP analyses of steady base flows that do not satisfy the steady Navier–Stokes equations, including timemean STBLIs (Touber & Sandham Reference Touber and Sandham2009; Pirozzoli et al. Reference Pirozzoli, Larsson, Nichols, Bernardini, Morgan and Lele2010; Waindim et al. Reference Waindim, Bhaumik and Gaitonde2016; Nichols et al. Reference Nichols, Larsson, Bernardini and Pirozzoli2017).
The timemean linear tendency is described in figure 8, visualized by the timemean streamwise velocity perturbation, for cases 11–13, in panels (a–c) respectively. We note that the mean perturbation field converges at a much slower rate than the mean base flow, since the range of the perturbation field extends several orders of magnitude. Thus, convergence of the mean is sensitive to the behaviour of large magnitude perturbations. As such, discrepancies between the nominal (a) and refined (b) grids are expected and observed, and they are attributable not only to differences in spatial refinement, but also to differences in initial STBLI realizations between different cases. Sensitivity to averaging duration can also occur, even with sufficiently long duration simulations which allow for accurate convergence of the timemean flow field; results on the nominal grid (a) are best converged in this sense, since case 11 has the longest simulation duration. This comparison also provides an indication of the level of uncertainty in the timemean perturbation quantities and sensitivity of these quantities to the base flow. Nonetheless, all cases, including the interaction with stronger separation (c, case 13), exhibit similar timemean linear tendencies, with equivalent physical significance.
The timemean behaviour of the absolute instability depicted in figure 8 is clarified by noting that, for all cases, in the timemean sense, the reflected shock position observed in the linearly constrained twin flow is upstream from its position in the base flow. This results in a negative timemean streamwise velocity perturbation in the region of the reflected shock as illustrated in figure 8(d); i.e. the mean linearized shock position is not coincident with the mean shock position, but is rather located slightly upstream. Notably, the unsteadiness of the base flow inhibits a sharp indication of the position of the linearized and baseflow shocks in the mean, since the difference in position between the shocks, at any given instant, is determined by the perturbation target magnitude, $\mathit{O}(10^{6})$ , which is certainly much less than the magnitude of unsteady shock displacement in the base flow, $\mathit{O}(1)$ ; however, the spatially diffuse negative timemean streamwise velocity perturbation in the region of the spatially diffuse timemean reflected shock reveals this behaviour. In conjunction with the timelocal analysis of § 5.3.6, we can furthermore conclude, regarding lowfrequency motion, that when the reflected shock is located near its mean position, the corresponding timelocal linear tendency promotes upstream shock motion, increasing the size of the separation region, consistent with the observed timemean linear tendency of the interaction. These findings are not inconsistent with those of Touber & Sandham (Reference Touber and Sandham2009), Pirozzoli et al. (Reference Pirozzoli, Larsson, Nichols, Bernardini, Morgan and Lele2010) and Nichols et al. (Reference Nichols, Larsson, Bernardini and Pirozzoli2017), who employ the timemean base flow to describe a linearly unstable nonoscillatory global mode. However, the current analysis shows that the tendency for upstream motion is definite; there is not an arbitrary nature to the sign (direction) of the absolute instability, it is not symmetric as suggested in previous modal meanflow analyses, and this observation has physical significance. Indeed, the current findings are consistent with the asymmetry of shock motion discussed by Piponniau et al. (Reference Piponniau, Dussauge, Debieve and Dupont2009), in which upstream motion occurs more rapidly than downstream motion. We discuss this in detail in § 5.3.5, followed by an analysis of how this absolute instability interacts with the nonlinear nature of the base flow in § 5.3.6.
We conclude this section by postulating that for any flow admitting a global absolute linear instability, the timemean linear tendency of the unsteady turbulent base flow should be qualitatively consistent with traditional stability analyses applied to the timemean turbulent flow. This postulate is confirmed in the current case of the STBLI, since the timemean linear tendency is consistent with the behaviour of the absolutely unstable, nonoscillatory, global mode of the timemean STBLI. Compared to stability analyses of the timemean turbulent flow, however, the dynamic linear response yields additional insight into the underlying mechanisms that sustain the unsteadiness, particularly in addressing timelocal perturbation dynamics, as discussed in § 5.3.6.
5.3.4 Effects of stochastic forcing on perturbation field
The linearly constrained perturbation field of the STBLI shows relatively small sensitivity to the specific nature of forcing, as anticipated given the absolute nature of the instability. Results with stochastic forcing are therefore similar to those shown earlier. This is consistent with the discussion by Chomaz (Reference Chomaz2005): in absolutely unstable flows, the internal growth of perturbations overshadows the response to external forcing. On the other hand, the linearly constrained perturbation field of the turbulent boundary layer, without an impinging shock, requires continuous stochastic forcing to sustain its convectively unstable perturbation field, consistent with the common understanding (Chomaz Reference Chomaz2005).
To demonstrate the selective nature of the STBLI perturbation field, the effects of forcing location on the rootmeansquare (r.m.s.) vertical velocity perturbation fields are displayed in figure 9. Apart from a local area of influence around each forcing location, the character of the perturbation field is very similar between the different cases. For instance in case 9, the forcing location provides a perturbation source upstream, but these perturbations are attenuated relatively quickly, and the major area of perturbation activity remains in the separation region and downstream boundary layer. Likewise, for case 7, the forcing effects are localized to the separated shear layer. For case 10, in which a perturbation source is added in the recirculation region, the results are quite similar to case 11, containing no forcing. This suggests that the selfsustaining perturbation field of case 11 is akin to a perturbation source originating in the recirculation region; i.e. the absolute instability is sustained by internal perturbation growth within the recirculation region. Additionally, it demonstrates that the STBLI is more sensitive to perturbations in the recirculation region than the upstream region; i.e. all attempts to stochastically force perturbations outside of the recirculation region are overshadowed by unforced perturbation growth in the recirculation region. These observations support the finding that the absolute instability is sustained by constructive feedback through the recirculation region and is not receptive to external forcing. The selfsustaining perturbation field, driven by the absolute instability, is the most natural perturbation field to analyse, since it removes any arbitrariness due to the presence of forcing, and is therefore the subject of much of the analysis below.
5.3.5 Spectral content of baseflow and perturbation fields
We now analyse the unsteady dynamics of the baseflow and perturbation fields using several different time series. These series are chosen to describe a variety of unsteady features: the shock footprint indicated by wall pressure fluctuations, indicators of the size of the separation bubble, as well as the nearfield shock location. Observations of the base flow are first summarized to provide a foundation for the subsequent discussion of the accompanying perturbation dynamics. Overall, the dynamics of the base flow agrees well with those described in the literature, providing a measure of baseflow validation necessary for the perturbation analysis. The discussion focuses on several prominent frequency bands, each described qualitatively with a corresponding physical process: lowfrequency content $(St_{L}\sim 0.03)$ , corresponding to separation bubble breathing and largescale shock motion, low–midfrequency content $(St_{L}\sim 0.1)$ , corresponding to bubble oscillations, high–midfrequency content $(St_{L}\sim 0.5)$ , corresponding to Kelvin–Helmholtz induced shear layer phenomena, including shedding and flapping of the reflected shock and highfrequency $(St_{L}\gtrsim 1)$ jitter, corresponding to the passage of finescale turbulent fluctuations.
In all cases, the time series are sampled at intervals of $20\unicode[STIX]{x1D6FF}_{0}/U_{\infty }$ , and the power spectral density (PSD) is obtained using Welch’s windowed overlapping segments estimator method, by splitting the signal into eight segments with 50 % overlap and windowing with the Hamming function. The PSD is premultiplied by the corresponding Strouhal number, based on the timemean separation length, to effectively represent the power distribution on a logarithmic scale. Further, the PSD is normalized, such that the integral of premultiplied power over all resolved frequencies is unity. In some cases, the PSD is filtered using Konno–Ohmachi (KO) smoothing, which ensures constant bandwidth smoothing on a logarithmic scale (Konno & Ohmachi Reference Konno and Ohmachi1998). The baseflow and perturbation fields were evaluated, in multiple simulations, for over five cycles of the lowfrequency $(St_{L}\sim 0.03)$ unsteadiness. For the former, adequacy was established by comparing with the known spectra from welldocumented results for this configuration (Dupont et al. Reference Dupont, Haddad and Debiève2006; Touber & Sandham Reference Touber and Sandham2009; Aubard et al. Reference Aubard, Gloerfelt and Robinet2013). With regard to the perturbations, these reach statistical stationarity quickly, and short duration simulations are observed to be adequate to describe their timelocal behaviour. For example, in principle, the bandisolation study (§ 5.3.6) only requires one cycle of the lowfrequency unsteadiness; results were nonetheless confirmed by comparing with the full data set. Therefore, the resolution of lowfrequency unsteadiness is reasonable to support the conclusions drawn in this work.
Notably, asymmetry is observed in the shock motion, whereby largedisplacement upstream motion occurs relatively quickly, followed by slower, less coherent, downstream motion, at intervals representative of the lowfrequency $(St_{L}\sim 0.03)$ cycle. This is evident in the baseflow nearfield reflected shock location signal (figure 10 a), which describes the streamwise location of the maximum value of velocity dotted with pressure gradient (constrained to the nearfield region of the reflected shock), at a height of $y=4\unicode[STIX]{x1D6FF}_{0}$ , in the spanwise centre of the domain. As noted earlier, this asymmetry is consistent with the experimental observations of Piponniau et al. (Reference Piponniau, Dussauge, Debieve and Dupont2009). Further analysis of the genesis of this asymmetry is provided in § 5.3.6.
Two additional signals are analysed to elucidate the dynamics of the separation region. A reverse flow massperspan signal, comprising the total fluid mass per span in the reversed flow region, is presented in figure 10(b), and a separation bubble massperspan signal, of the total fluid mass per span satisfying $u<0.15u_{\infty }$ , is presented in figure 10(c). The reverse flow massperspan signal (figure 10 b) is dominated by $St_{L}\sim 0.1$ content, corresponding to bubble oscillations, with a minor side lobe at $St_{L}\sim 0.03$ , corresponding to lowfrequency bubble breathing. On the other hand, the separation bubble massperspan signal: $u<0.15u_{\infty }$ (figure 10 c) shows a smaller contribution from $St_{L}\sim 0.1$ content, and a larger contribution from $St_{L}\sim 0.03$ content, resulting in a better indication of shock displacement. This trend holds across other thresholds defining the ‘turbulent separation bubble’: larger thresholds on streamwise velocity result in a better indication of lowfrequency $(St_{L}\sim 0.03)$ largescale shock motion, while smaller thresholds result in a better indication of low–midfrequency $(St_{L}\sim 0.1)$ bubble oscillation dynamics. As § 5.3.6 elaborates, the dynamics in these separate frequency bands can be described relatively independently. Priebe & Martín (Reference Priebe and Martín2012) similarly observe significant $St_{L}\sim 0.1$ content corresponding to streamwise oscillation of the lowpass filtered separation point in the DNS of a Mach 2.9 compression ramp STBLI. In fact, $St_{L}\sim 0.1$ content is found in the separation region power spectra of many studies (Dupont et al. Reference Dupont, Haddad and Debiève2006; Touber & Sandham Reference Touber and Sandham2009; Aubard et al. Reference Aubard, Gloerfelt and Robinet2013), but it is often not discussed in detail.
We also observe a significant lagged correlation between the reflected shock motion and the dynamics of the separation bubble. The maximum correlation coefficient between the shock location and the reverse flow massperspan signals is 0.38 with a lag of $18\unicode[STIX]{x1D70F}$ . A stronger instantaneous correlation of 0.53 is found with a lag of $15\unicode[STIX]{x1D70F}$ , when correlating the shock location and separation bubble massperspan: $u<0.15u_{\infty }$ signals. These correlations of shock motion with the size of the separation, in which the shock motion lags in response to changes in separation, are consistent with previous observations by Morgan et al. (Reference Morgan, Duraisamy, Nguyen, Kawai and Lele2013), confirming that a causal relationship could exist, whereby largescale shock motion responds primarily to growth or contraction of the separation bubble. This causal relationship (sign of the lag) is independent of the height at which the correlation with the reflected shock is taken, though the magnitude of the correlation and lag vary.
In figure 11, the base flow and perturbation wall pressure PSDs are compared at streamwise locations around the interaction region. The PSD at each streamwise location is normalized, such that the integral of premultiplied power over all resolvable frequencies is unity. The baseflow wall pressure dynamics agree well with previous studies (Dupont et al. Reference Dupont, Haddad and Debiève2006; Touber & Sandham Reference Touber and Sandham2009; Agostini et al. Reference Agostini, Larchevêque, Dupont, Debiève and Dussauge2012; Aubard et al. Reference Aubard, Gloerfelt and Robinet2013). Lowfrequency content $(St_{L}\sim 0.03)$ , corresponding to shock motion and bubble breathing, is prominent near the timemean separation and reattachment locations; normalizing the PSD independently at each streamwise location facilitates the identification of prominent lowfrequency dynamics at reattachment; these fluctuations are less intense than the lowfrequency fluctuations at separation, but they are intense compared to the pressure fluctuations across the full spectrum at reattachment. Low–midfrequency content $(St_{L}\sim 0.1)$ , corresponding to bubble oscillations, is present near the separation and reattachment points, as well as downstream of the interaction region. High–midfrequency content $(St_{L}\sim 0.5)$ , corresponding to Kelvin–Helmholtz shedding, is present downstream of the interaction region. Highfrequency content $(St_{L}\gtrsim 1)$ , typical of finescale features in a turbulent boundary layer, is prominent upstream of the interaction region; it persists throughout the interaction, reducing in frequency as its convective speed slows and its outer scale increases, eventually being overtaken in amplitude by a dominant contribution from $(St_{L}\sim 0.5)$ shed structures downstream of the interaction.
Necessarily, the unsteady dynamics in each of these frequency bands will not collapse with $St_{L}=fL/U_{\infty }$ for all flow conditions. Rather, $St_{L}$ is chosen to describe the unsteadiness since the prominent lowfrequency unsteadiness collapses with this Strouhal number for a variety of flow conditions (Dussauge et al. Reference Dussauge, Dupont and Debiève2006), providing a unified scale on which to compare all of the described bands with respect to the lowfrequency unsteadiness for these particular flow conditions. Indeed, the highfrequency band associated with finescale turbulent fluctuations of the incoming boundary layer must collapse with $St_{\unicode[STIX]{x1D6FF}}=f\unicode[STIX]{x1D6FF}_{0}/U_{\infty }=1$ , independent of the STBLI separation length. It is still to be determined for which Strouhal numbers the two described midfrequency bands might collapse, since these frequency bands have not been characterized in as extensive detail in many studies.
The wall pressure PSD for the perturbation field is calculated in the same manner as for the base flow. Downstream of the interaction, high frequencies typical of a turbulent boundary layer $(St_{L}\gtrsim 1)$ are identified, consistent with the observation that the perturbation field highlights regions of highpotential linear growth of the baseflow turbulent structures. Midfrequency $(St_{L}\sim 0.060.1)$ content is found near the timemean separation point, corresponding to perturbations of the reflected shock. This suggests that changes in the linear tendency of the reflected shock are more significantly influenced by midfrequency mechanisms than lowfrequency mechanisms. The wall pressure perturbation PSD upstream of the timemean separation point requires special care in its interpretation, since the perturbation field does not always penetrate much upstream of the timemean separation point. For instance, if a particular upstream lowfrequency motion of the reflected shock does not reach its greatest upstream extent, the perturbation field will also not be carried farther upstream. The effect of this is to skew the PSD toward higher frequencies. Nonetheless, this analysis illustrates that while the perturbation field is a strong function of the base flow, there are significant differences in the spectral properties of the baseflow and perturbation fields.
To demonstrate the linearity of forcing and attenuation on perturbation field spectral content, wall pressure spectra are presented in figure 12, for cases 6–8, in which forcing and attenuation are applied with different magnitudes, separated by one order of magnitude, as described in table 4. In these spectra, the PSD is normalized by $(\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P_{L\,target}})^{2}$ , under which the spectra should collapse if linearity is maintained. KO smoothing has been applied to reduce noise at higher frequencies. At the timemean separation point (a), cases 7 and 8 collapse, while case 6 undershoots the peak. The collapse of cases 7 and 8 is anticipated; however, the deviation of case 6 appears to demonstrate nonlinearity in the scaling of forcing and attenuation target magnitudes between cases. This apparent nonlinearity can result from differences in the specific realization of stochastic forcing, which is unique to each case; it is not indicative that higherorder terms in the expansion of (3.11) are nonnegligible, since the results § 5.3.2 demonstrate that the scaling of attenuation alone is linear. Similarly, the results of Unnikrishnan & Gaitonde (Reference Unnikrishnan and Gaitonde2016) demonstrate that the method is linear with respect to scaling of deterministic forcing. Regardless, at the downstream location (b), all cases collapse indicating linearity. The results presented herein regarding the selfsustaining linearly constrained perturbation field are obtained with the attenuation factor of case 7, reaffirming that the perturbations remain linear with respect to the base flow.
The attenuation factor is related to the enforcement of linearity and thus provides insight into the dynamics of the interaction, since it distinguishes instances in time where perturbation growth is more prominent from those where it is less so. Characteristics of the variable attenuation factor for the selfsustaining perturbation field (case 11) are presented in figure 13. The attenuationfactor PSD is characterized by two peaks at approximately $St_{L}\sim 1$ , corresponding to the passage of convecting finescale turbulent structures, and $St_{L}\sim 200$ , corresponding to the changing location of the maximum norm of the perturbation field $(\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P})$ with the finescale features of convecting structures; recall, the perturbation field evolves to highlight turbulent features with high potential for linear growth. Scrutiny of animations of the evolving perturbation and baseflow fields indicates that the maximum norm of the perturbation field $(\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P})$ is often associated with a hairpinlike structure downstream of the interaction region. Since the attenuation factor, $\unicode[STIX]{x1D6FC}$ , is directly dependent on $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P}$ , it varies with the coherence time of these turbulent structures, explaining the lowfrequency peak. Additionally, the band of frequencies comprising this peak indicates that highfrequency $(St_{L}\sim 1)$ turbulent content and high–midfrequency $(St_{L}\sim 0.1{}0.5)$ content have a much stronger tendency for linear growth in the baseflow turbulent environment than lowfrequency $(St_{L}\sim 0.03)$ content. It follows that any mechanism for instability acting over lowfrequency scales must have a much smaller growth rate than its higherfrequency $(St_{L}\sim 0.1{}2)$ counterparts. The highfrequency $(St_{L}\sim 200)$ peak corresponds to variation of $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P}$ with shorttime evolution of finescale features on a passing coherent structure. Since the turbulent structure evolves throughout its time of coherence, the spatial location at which $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P}$ is found will vary with the evolution of that turbulent structure as the location for maximum potential growth changes. The attenuation factor is smoothly resolved by the chosen time step, with the highfrequency variation of $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}}^{P}$ being resolved by $\mathit{O}(10)$ time steps. Note that the discrete attenuation factor $(1\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}\unicode[STIX]{x1D70F})$ exceeds unity in some instances, leading to brief amplification of the perturbation field. On average, however, the discrete attenuation factor takes a value of 0.860, with a standard deviation of 0.084. It follows that, on average, the absolute instability of the STBLI corresponds to a Lyapunov exponent of approximately $140U_{\infty }/\unicode[STIX]{x1D6FF}_{0}$ , significantly larger than the Lyapunov exponent of $10.3U_{\infty }/\unicode[STIX]{x1D6FF}_{0}$ identified earlier for the shockfree turbulent boundary layer; i.e. the STBLI absolute instability results in an exponential rate of perturbation amplification that is an order of magnitude larger than the shockfree, convectively unstable, turbulent boundary layer.
5.3.6 Band isolation of baseflow and perturbation fields
In this section, we isolate the flow characteristics associated with the prominent frequency bands observed in the spectral analysis of § 5.3.5. For this purpose, a series of lowpass finiteimpulseresponse (FIR) windowbased filters are applied. Since the lowfrequency dynamics are broadband, a sharp cutoff is neither necessary nor desirable. They are therefore implemented through a time domain moving average filter, to optimize smoothing characteristics in the time domain, albeit with a wider transition band. The first lowpass cutoff is chosen at $St_{L}=0.06$ , to isolate lowfrequency unsteadiness of bubble breathing around $St_{L}\sim 0.03$ . The lowpass content is then subtracted from the instantaneous signal, and a second lowpass filter is applied with a cutoff at $St_{L}=0.2$ , to isolate the low–midfrequency bubble oscillations around $St_{L}\sim 0.1$ . The first two modes are then subtracted from the instantaneous signal, and a third filter with lowpass cutoff frequency of $St_{L}=0.7$ is applied, to isolate high–midfrequency Kelvin–Helmholtz shedding around $St_{L}\sim 0.5$ . The residual contains higher frequencies, typical of convecting finescale turbulent structures, discussed earlier in the context of figure 13. This procedure results in a set of bandpass filtered modes that form a complete basis for the timeresolved flow.
This technique allows for band isolation of the flow similar to previous analyses of STBLIs. Several efforts (Pirozzoli et al. Reference Pirozzoli, Larsson, Nichols, Bernardini, Morgan and Lele2010; Aubard et al. Reference Aubard, Gloerfelt and Robinet2013; Priebe et al. Reference Priebe, Tu, Rowley and Martín2016; Nichols et al. Reference Nichols, Larsson, Bernardini and Pirozzoli2017) have accomplished this by applying a Fourier decomposition, or DMD, to fluctuations of a stationary flow; the latter approach is equivalent to the former under these conditions (Chen et al. Reference Chen, Tu and Rowley2012). Others (Pirozzoli et al. Reference Pirozzoli, Larsson, Nichols, Bernardini, Morgan and Lele2010; Priebe & Martín Reference Priebe and Martín2012) have similarly employed filtering techniques for this purpose. As anticipated, bandisolated modes of the STBLI agree qualitatively well, regardless of the technique used to generate them. For this application, in which the frequencies of interest are broadband (figure 11 a), we employ the filtering approach, from which a single mode reproduces the timelocal dynamics more faithfully. The bandisolated modes of the base flow agree well with the trends observed in these previous studies; the lowfrequency $(St_{L}\sim 0.03)$ and high–midfrequency $(St_{L}\sim 0.5)$ modes have been well documented and compare favourably. We emphasize that the primary insight provided in this section is not to validate the observed bandisolated behaviour of the base flow, but rather to band isolate the base flow and perturbations simultaneously, drawing insight from how the bandisolated perturbations are correlated with the base flow at certain phases of each bandisolated cycle. This analysis provides unique information about linear tendencies, which was not available to the aforementioned efforts that analysed the bandisolated base flow.
Figures 14, 16 and 17 display a representative cycle of the bandpass filtered results corresponding to each of these frequency bands. In each, the left column shows the bandpass filtered streamwise velocity baseflow component, with the timemean streamwise velocity added to aid in visualization, while the right column shows the bandpass filtered streamwise velocity perturbation component. Each cycle is described by snapshots taken at evenly spaced intervals. Contour levels are consistent through each column in all figures, and the spatial coordinates ( $x$ – $y$ plane) are normalized by the timemean separation length with the streamwise origin located at the timemean separation point.
Just as the timemean linearly constrained perturbation quantities represent the timemean linear tendency of the base flow, the filtered (shorttimemean) linearly constrained perturbation quantities represent the timelocal linear tendency of the base flow in the corresponding frequency band. This ability to probe the timelocal linear tendency of the base flow is one of the primary advantages of the LCSLES method over perturbation analyses of the timemean flow. Provided the filter cutoff frequency is moderately lower than any turbulent noise (i.e. ensuring the flow evolves long enough for ergodic averaging to smooth out turbulent noise), the resulting timelocal linear tendencies can provide relatively clean indications of the timelocal character of flow stability. For instance, one may examine the timelocal linear tendency associated with the baseflow reflected shock at an upstream position and infer, at that instant, the direction in which linear mechanisms are driving the baseflow reflected shock. In higherfrequency bands, the timelocal linear tendencies are somewhat more contaminated by highfrequency turbulence due to a limited duration for ergodic averaging, as well as flow history effects (i.e. due to the hyperbolic – wavelike – nature of the perturbations in time, the past history of the perturbations becomes more significant, and the observed linear tendencies will lag behind the causal mechanisms in the base flow at higher frequencies), which must be factored into the analysis.
Breathing of separation bubble $(St_{L}\sim 0.03)$ . Snapshots from the lowpass filtered $(St_{L}<0.06)$ flow field are shown in figure 14. The results are discussed by considering four phases of the breathing cycle at frequency $St_{L}\sim 0.03$ .
 $\mathit{Phase}=0:$

The bubble volume is at a neutral extent and in the depletion phase of the lowfrequency cycle. Negative $u^{P}$ perturbations along the separated shear layer and in the region of the reflected shock indicate the linear tendency for retardation of the flow, upstream motion of the shear layer/shock and bubble growth. Nonetheless, the bubble continues to contract due to nonlinear massdepletion mechanisms discussed below.
 $\mathit{Phase}=\unicode[STIX]{x03C0}/2:$

The bubble volume is at its smallest extent. Like phase 0, negative $u^{P}$ perturbations indicate the linear tendency for bubble growth, and the magnitude of (negative) $u^{P}$ in this phase is larger than for phase 0, suggesting a greater potential for growth than before. The bubble now begins to grow.
 $\mathit{Phase}=\unicode[STIX]{x03C0}:$

The bubble volume is at a neutral extent and in the growth phase of the lowfrequency cycle. Negative $u^{P}$ perturbations similarly indicate the linear tendency for growth, and the bubble continues to grow.
 $\mathit{Phase}=3\unicode[STIX]{x03C0}/2:$

The bubble volume is at its largest extent. Now positive $u^{P}$ perturbations are prominent along the separated shear layer and in the region of the reflected shock, indicating the linear tendency for acceleration of the flow, downstream motion of the shear layer/shock, and bubble depletion. The bubble begins to contract.
The linearperturbation state (figure 14 b,d,f,h) is consistent with the observed motion of the separation bubble and reflected shock in the base flow (figure 14 a,c,e,g) at phases $\unicode[STIX]{x03C0}/2$ , $\unicode[STIX]{x03C0}$ , and $3\unicode[STIX]{x03C0}/2$ . This suggests that at these times, the bubble/shock system is responding to the timelocal linear tendency of the base flow. However, at phase 0, the perturbation state is contrary to the motion of the bubble, and the bubble continues to contract, while the perturbations indicate that it should grow. This suggests that a nonlinear mechanism is responsible for bubble contraction at this time. This nonlinear depletion process could be due to the mass entrainment mechanism proposed by Piponniau et al. (Reference Piponniau, Dussauge, Debieve and Dupont2009), in which fluid shedding from the recirculation region accounts for the nonlinear dynamics driving bubble contraction. We further note that while this lowfrequency bubble/shock motion occurs at approximately periodic intervals throughout the simulation, it is not generally distributed uniformly throughout each cycle. In fact, the observed motion is similar to the experimental observations of Piponniau et al. (Reference Piponniau, Dussauge, Debieve and Dupont2009), in that the bubble remains at a small or moderate size through most of the cycle, followed by rapid growth, then a moderate rate of bubble depletion accompanied by the corresponding shock motion. The cycle asymmetry is not necessarily evident in the lowpass phase diagrams, as the sudden growth is smeared from phase $3\unicode[STIX]{x03C0}/2$ into phase $\unicode[STIX]{x03C0}$ , resulting in a larger bubble at phase $\unicode[STIX]{x03C0}$ than would typically be observed in the instantaneous flow. However, the lineartendency trends represented in this cycle are informative in describing the origin of cycle asymmetry.
In fact, the observed behaviour is consistent with interplay between the baseflow linear tendency and massdepletion (nonlinear) responses of the bubble. In the initial phases of the cycle, the linear tendency of the base flow promotes bubble growth, while mass depletion promotes contraction. As such, the bubble is suspended in a state of near equilibrium for a time, as mass depletion gradually reduces the bubble size through phases $\unicode[STIX]{x03C0}/2$ – $\unicode[STIX]{x03C0}$ . At this point, mass depletion from the bubble nears completion, so nonlinear forcing promoting bubble contraction diminishes, and the baseflow linear tendency promoting growth overcomes the contraction effect of mass depletion. The bubble then grows rapidly to phase $3\unicode[STIX]{x03C0}/2$ ; i.e. the imbalance, resulting from the reduction in downstream nonlinear forcing relative to upstream linear tendency, allows for the sudden and intense upstream mass flux in the reverse direction, during which the bubble grows to the state observed at phase $3\unicode[STIX]{x03C0}/2$ . However, the growth is abruptly slowed once the character of the baseflow linear tendency changes as shown at phase $3\unicode[STIX]{x03C0}/2$ , at which point the bubble is so large that the baseflow linear tendency now promotes contraction. The upstream motion of the reflected shock is reversed, and the bubble quickly relaxes to a more moderate state. The cycle then resets, and the bubble is suspended in a state of near equilibrium due to the competing lineartendency and massdepletion mechanisms. While the essential components of this mechanism have been proposed previously (Plotkin Reference Plotkin1975; Piponniau et al. Reference Piponniau, Dussauge, Debieve and Dupont2009), the current dynamic linear response, obtained using the LCSLES technique, suggests a synthesis of these concepts, providing the first demonstration of causal dynamics embedded in the timelocal linear tendency of the baseflow lowfrequency band. It also facilitates the description of a possible mechanism capable of explaining asymmetry in the shock motion cycle, of which the key component is found in the (phase 0) competing mechanisms of baseflow linear tendency and (nonlinear) mass depletion. Notably, the above description may be extended, in the case of an STBLI configuration with different baseflow properties, to explain the opposite asymmetrical bias; e.g. if the baseflow timemean linear tendency described a downstream tendency for shock motion, yet the interaction retained linear restoring tendencies at extreme shock displacements, then the above explanation may be applied, with inverted directional bias, leading to rapid downstream shock motion and contraction of the separation bubble, followed by moderate upstream shock motion and bubble growth; however, such behaviour is not observed in any of the cases investigated in this work.
This analysis also suggests that the STBLI permits a timemean linearly stable configuration, in which the reflected shock is slightly upstream from the timemean location. At the extremes of the lowfrequency cycle (phases $\unicode[STIX]{x03C0}/2$ and $3\unicode[STIX]{x03C0}/2$ ), the bubble and shock motion are consistent with a linear tendency for restoration to a more moderate configuration; previously it was demonstrated that the timemean base flow is linearly unstable, with a timemean linear tendency for upstream reflected shock motion, similar to the linearly unstable nonoscillatory mode found in the IVP calculations of Touber & Sandham (Reference Touber and Sandham2009) and the IVP/EVP calculations of Pirozzoli et al. (Reference Pirozzoli, Larsson, Nichols, Bernardini, Morgan and Lele2010) and Nichols et al. (Reference Nichols, Larsson, Bernardini and Pirozzoli2017). Hence, the linearly stable configuration would consist of a separation bubble slightly larger than the timemean size and a reflected shock slightly upstream from the timemean location, as indicated in figure 15. Lowfrequency unsteadiness could then be described in terms of nonlinear mechanisms forcing the reflected shock, which at extreme displacements experiences a linear tendency for restoration to more moderate displacements, as is the premise of the ODE model proposed by Plotkin (Reference Plotkin1975) and derived by Touber & Sandham (Reference Touber and Sandham2011). However, we find that the ‘timemean linearly stable shock position’ and the ‘timemean shock position’ do not coincide, contrary to the assumption of this ODE model, a correction for which could lead to an improved model. We suggest implementing such a correction in the form of a damped and driven harmonic oscillator model, as opposed to the Langevin model that is currently employed; necessarily, to reproduce the shock motion asymmetry with such a model, the stochastic forcing that drives the oscillator should have a functional dependence on the position and velocity of the shock, reflecting contributions from internal mechanisms; e.g. the downstream forcing, resulting from bubble mass depletion, must diminish as the bubble contracts. It follows that, in addition to common techniques which reduce the size of the separation, concepts for control that reduce the discrepancy between the ‘timemean’ and ‘timemean linearly stable’ shock positions may prove effective at muting lowfrequency unsteadiness.
Oscillation of separation bubble $(St_{L}\sim 0.1)$ . Figure 16 shows snapshots of the bandpass filtered $(0.06<St_{L}<0.2)$ flow field highlighting oscillation of the separation bubble at frequency $St_{L}\sim 0.1$ . In this band, six phases are employed to describe the cycle.
 $\mathit{Phase}=0:$

The bubble (white oval) is downstream (position $(0.7,0.02)$ ) moving upstream. The separated shear layer linearly tends to retard upstream motion (indicated by positive $u^{P}$ ), but along the wall (position $(0.5,0.02)$ ), there is a linear tendency for upstream acceleration of bubble motion (indicated by negative $u^{P}$ upstream of the current bubble location).
 $\mathit{Phase}=\unicode[STIX]{x03C0}/3:$

The bubble is midstream (position $(0.5,0.02)$ ) moving upstream. The shear layer linearly tends to retard upstream motion, but along the wall (position $(0.5,0.02)$ ), there is a linear tendency for acceleration of upstream bubble motion (more extreme than during phase 0).
 $\mathit{Phase}=2\unicode[STIX]{x03C0}/3:$

The bubble is upstream (position $(0.2,0.02)$ ), and accelerated by the linear tendency along the wall (position $(0.2,0.02)$ ), but decelerated by the linear tendency away from the wall (position $(0.2,0.1)$ ). The bubble splits: low velocity fluid (pink oval) moves away from wall and is shed downstream, accelerated by the shear layer linear tendency (position $(0.2,0.1)$ ). The shed fluid continues moving downstream.
 $\mathit{Phase}=\unicode[STIX]{x03C0}:$

The bubble moves relatively quickly downstream with the shed fluid (position $(0.7,0.02)$ ). This motion is retarded by the shear layer linear tendency.
 $\mathit{Phase}=4\unicode[STIX]{x03C0}/3:$

The bubble (position $(0.9,0.02)$ ) and shed fluid continue to move downstream. This downstream motion continues to be retarded by the shear layer linear tendency.
 $\mathit{Phase}=5\unicode[STIX]{x03C0}/3:$

The bubble reaches its farthest downstream extent (position $(1.0,0.02)$ ). The bubble and shed fluid are now accelerated downstream, out of interaction region, in accordance with the observed linear tendency. Separation size achieves a minimum. A new nearwall kernel with linear tendency promoting bubble growth and upstream motion appears (position $(0.4,0.02)$ ).
This cycle demonstrates the oscillatory dynamics of the separation bubble at $St_{L}\sim 0.1$ . At phase 0, a small negative $u^{P}$ kernel appears near the midpoint of the interaction region. The bubble emerges near this position and has a linear tendency to be driven upstream, along the wall, in the initial phases $0{}\unicode[STIX]{x03C0}$ , as demonstrated by the local region of negative $u^{P}$ upstream of the current bubble position. During this upstream bubble motion, the linear tendency of the shear layer opposes upstream motion, except in the nearwall, nearbubble region, as demonstrated by a positive region of $u^{P}$ . After traversing upstream to near the timemean separation point, the bubble splits, and some low velocity fluid moves away from the wall entering the offwall region of faster moving flow (phase $2\unicode[STIX]{x03C0}/3$ ). The bubble then quickly moves downstream accompanying the shed fluid (phases $\unicode[STIX]{x03C0}{}4\unicode[STIX]{x03C0}/3$ ). The shear layer linearly tends to oppose this downstream motion, as indicated by a broad region of negative $u^{P}$ . Finally, after traversing downstream to near the timemean reattachment point, the shed fluid now linearly tends to be accelerated out of the interaction region (indicated by a region of positive $u^{P}$ ), and the bubble contracts. A new kernel of negative $u^{P}$ appears near the interaction region midpoint, and the cycle repeats.
This behaviour is qualitatively similar to several weakly damped oscillatory modes described by Pirozzoli et al. (Reference Pirozzoli, Larsson, Nichols, Bernardini, Morgan and Lele2010) and Nichols et al. (Reference Nichols, Larsson, Bernardini and Pirozzoli2017) in their BiGlobal stability analyses of timemean STBLI flow fields. Specifically, Nichols et al. (Reference Nichols, Larsson, Bernardini and Pirozzoli2017) describes several weakly damped oscillatory modes $(0.06<St_{L}<0.13)$ for a Mach 2.28 STBLI, with an impinging shock generated by a $10.7^{\circ }$ farfield wedge. These modes depict periodic oscillations of the separation bubble and reflected shock, similar to our observations of the bandpass filtered $(St_{L}\sim 0.1)$ flow fields. The ‘type II’ modes of Nichols et al. (Reference Nichols, Larsson, Bernardini and Pirozzoli2017), in which the ‘reflected shock plays a more passive role’, are particularly similar to the phases of the separation bubble described in figure 16. Additionally, $(St_{L}\sim 0.1)$ content is found in the reflected shock signal (figure 10 d), consistent with the ‘type I’ modes, of the BiGlobal analysis, which primarily depict periodic reflected shock motion. However, we suggest a new interpretation of these observations, in conjunction with our recent observations of the dynamic linear response of the timeresolved STBLI: associating the zerofrequency (nonoscillatory) unstable mode with lowfrequency $(St_{L}\sim 0.03)$ unsteadiness, as discussed in the context of § 5.3.3 and figure 14, and associating the weakly damped oscillatory modes with low–midfrequency $(St_{L}\sim 0.1)$ unsteadiness, as discussed in the context of figure 16. Nonetheless, qualitative agreement between the observed oscillatory behaviour in the LES and BiGlobal stability analyses confirm this mechanism to be an important aspect of low–midfrequency $(St_{L}\sim 0.1)$ STBLI dynamics.
The bandisolated dynamics of the base flow and perturbations are consistent with the action of a massdepletion mechanism, such as in the conceptual model of Piponniau et al. (Reference Piponniau, Dussauge, Debieve and Dupont2009); the competition of the accompanying nonlinear downstream forcing and upstream linear tendency are necessary to describe the lowfrequency dynamics. However, like Morgan et al. (Reference Morgan, Duraisamy, Nguyen, Kawai and Lele2013) we find that the ‘turbulent separation bubble’ (§ 5.3.5) responds more significantly at higher frequencies $(St_{L}\sim 0.1)$ than predicted by the model $(St_{L}\sim 0.03)$ . Higherfrequency $(St_{L}\sim 0.1)$ content is also prominent in the dynamics of subsonic separation unsteadiness (Kiya & Sasaki Reference Kiya and Sasaki1985), and it is not entirely clear that the difference in mixing layer spreading rate, between the low and highspeed separated flows, fully accounts for the lower frequency $(St_{L}\sim 0.03)$ content in the STBLI, without additional considerations of shockinduced effects. Indeed, Agostini, Larchevêque & Dupont (Reference Agostini, Larchevêque and Dupont2015) discuss the significance of the recirculation bubble, with effects unique to this STBLI configuration, in governing both low and midfrequency shock dynamics for moderately to massively separated interactions, drawing comparison with aspects of incompressible separated flows.
Furthermore, the prominence of this frequency band in our reverse flow and separation bubble spectra suggests that the accompanying dynamics (figure 16) may be the dominant mechanism for mass transport from the recirculation region to the boundary layer. However, this does not eliminate the possibility for entrainment associated with Kelvin–Helmholtz $(St_{L}\sim 0.5)$ shedding to contribute to mass transport from the recirculation region; ultimately, all of these processes are nonlinearly coupled. Notably, this oscillatory mechanism may also be responsible for the ‘lowfrequency’ unsteadiness observed by Kiya & Sasaki (Reference Kiya and Sasaki1985). As discussed previously, a potential mechanism for lowfrequency $(St_{L}\sim 0.03)$ unsteadiness in the STBLI relies on interplay between the linear tendency of the STBLI base flow and nonlinear forcing associated with this oscillatory $(St_{L}\sim 0.1)$ massdepletion mechanism. In the case of Kiya & Sasaki (Reference Kiya and Sasaki1985), there is no shock and thus no interplay to produce unsteadiness at $St_{L}\sim 0.03$ , in which case oscillatory bubble mass changes near the frequency $St_{L}\sim 0.1$ are the lowestfrequency unsteady features present in that flow. That is, the linear component (figure 14) of the mechanism proposed to describe STBLI unsteadiness at $St_{L}\sim 0.03$ only pertains to flows with a similar timemean linear tendency to the STBLI, and may not pertain to incompressible separated flows, since the separation is not shock induced. For these cases, the linear oscillatory mechanism (figure 16) proposed to describe STBLI unsteadiness at $St_{L}\sim 0.1$ may be the primary mechanism of ‘lowfrequency’ unsteadiness.
Kelvin–Helmholtz (K–H) shedding $(St_{L}\sim 0.5)$ . Figure 17 shows snapshots from the bandpass filtered $(0.2<St_{L}<0.7)$ flow field highlighting K–H phenomena at frequency $St_{L}\sim 0.5$ . The base flow clearly demonstrates wavelike motion through the recirculation region, which corresponds to flapping of the reflected shock foot. K–H waves in this frequency band contribute to the origin of shed structures, leading to dominant midfrequency $(St_{L}\sim 0.5)$ content in the postinteraction wall pressure spectra of figure 11(a). We note, however, as shown in the reverse flow and separation bubble spectra of figure 10(d), that the mass of fluid in the recirculation region does not respond significantly at this frequency, and instead simply rides the waves. This suggests that while K–H $(St_{L}\sim 0.5)$ waves originate in the separation region and account for the increase in convective content around $St_{L}\sim 0.5$ downstream of the interaction, it is insufficient to describe the midfrequency separation bubble dynamics without noting both the oscillatory behaviour in the $St_{L}\sim 0.1$ band, as well as wavelike behaviour in the $St_{L}\sim 0.5$ band.
Since this band approaches higher frequencies $(St_{L}\sim 1)$ associated with convecting finescale turbulence, we observe more significant contamination from turbulence and flow history effects in the corresponding bandpass filtered perturbation field than in previous lowerfrequency bands. Still, periodic changes in the perturbation field character near the timemean separation point and along the separated shear layer coincide with the genesis frequency and convection of K–H waves, suggesting that these may also be driven by a periodic linear mechanism acting throughout the interaction; the separation bubbles are approximately driven along the major gradients of high to low (red to blue) streamwise velocity perturbation, consistent with the observed linear tendency. The clarity of the bandpass filtered perturbation field in this higherfrequency band could be improved through a conditional analysis, but is not considered in the current work.
This sequential lowpass filtering analysis thus aids in understanding the dynamic relation between the baseflow and linearly constrained perturbation fields in different frequency bands. In each frequency band, and at specific phases, it allows for the identification of whether the STBLI is responding consistently with the timelocal linear mechanisms identified in the evolving flow, or responding contrary to the identified linear mechanisms, suggesting the increased importance of nonlinear mechanisms. We anticipate that these linear mechanisms are an essential source of unsteadiness in all similar STBLIs, including those that are influenced by lowfrequency coherent fluctuations of the incoming flow, though such coherent fluctuations are not observed in the present study; the relative significance of internal to external mechanisms as sources of unsteadiness remain flowspecific questions. Although this analysis does not demonstrate a statistically robust measure of the dynamical behaviour in each frequency band, it does identify clear trends associated with each frequency band that are qualitatively similar throughout many realizations. As such, it successfully illustrates several key unsteady broadband features present in the STBLI and aids in identifying potential linear mechanisms which may sustain these features in the unsteady turbulent environment.
6 Conclusions
A Mach 2.3 shock/turbulentboundarylayer interaction is studied through small amplitude (linear) forcing of the unsteady separated flow. The response to this forcing is obtained by solving a constrained linearization about the timeevolving turbulent flow. For this purpose, a synchronized pair of largeeddy simulations is employed, in which the two simulations are advanced synchronously in time, with the second simulation constrained linearly to the first. The two solutions are processed to yield the evolution of linear perturbations of the fully turbulent STBLI. Ensemble statistics of the linear perturbations constitute the dynamic linear response of the flow, analysis of which allows for the inference of both the timemean as well as the timelocal linear tendencies of the unsteady turbulent flow. In limiting cases, these tendencies are analogous to and qualitatively consistent with traditional stability analyses applied to the steady timemean base flow.
The dynamic linear response demonstrates that the STBLI fosters a global absolute linear instability corresponding to a linear tendency for upstream shock motion and enlargement of the separation bubble. The absolute instability is maintained through constructive feedback of perturbations through the recirculation; it is selfsustaining and insensitive to external forcing. In contrast, this phenomenon is not observed in similarly perturbed turbulent boundary layers, in which linearly constrained perturbation growth is convective in nature and sustained only through continuous external forcing. The selfsustaining perturbation field of the STBLI is similar in character to the perturbation field forced at a location inside the timemean recirculation region, indicating that perturbations originating in the recirculation region have a more significant influence on STBLI dynamics than those originating upstream. The perturbation field remains anchored to the interaction region and exhibits a timemean linear tendency for upstream motion of the reflected shock. This nonzero timemean linear tendency quantifies the nature of the global absolute linear instability, and its bias for upstream motion is inherently coupled to observed asymmetry in the lowfrequency shock motion cycle, in which upstream shock motion occurs rapidly followed by slower downstream motion.
Significant dynamic features of the STBLI are identified and delineated by key frequency bands: $St_{L}\sim 0.5$ corresponding to high–midfrequency Kelvin–Helmholtz shedding along the separated shear layer, $St_{L}\sim 0.1$ corresponding to low–midfrequency oscillations of the separation bubble, and $St_{L}\sim 0.03$ corresponding to lowfrequency largescale bubble breathing and shock motion. Informed by observation of these frequencies in the STBLI, a bandpass filtering decomposition is applied to the STBLI baseflow and perturbation fields to isolate these broadband features and identify how the timelocal linear tendencies of the STBLI contribute to midfrequency and lowfrequency dynamics.
The lowfrequency $(St_{L}\sim 0.03)$ breathing of the separation bubble and corresponding shock motion can be explained in the context of competing linear and nonlinear mechanisms; linear mechanisms are identified from the dynamic linear response, and nonlinear mechanisms include mass depletion from the recirculation region through shedding during bubble contraction. This proposed explanation also accounts for the observed asymmetry in the lowfrequency cycle. Furthermore, this analysis confirms the presence of linear restoring tendencies in the STBLI; when the reflected shock is located at an extreme upstream or downstream displacement from its timemean position, a linear tendency for restoration to a more moderate position is observed. However, this moderate position (a linearly stable position in the timemean sense) does not correspond with the timemean shock position, but rather is located upstream from the timemean shock position. This observation that the ‘timemean linearly stable position’ and the ‘timemean position’ do not coincide supports the assertion that the lowfrequency unsteadiness is an intrinsic property of the system, which could potentially be mitigated through control by altering the interaction in a way such that the two positions are more closely coincident.
The separation bubble is found to oscillate at midfrequency $(St_{L}\sim 0.1)$ in accordance with the linear tendency of the STBLI, and the mass budget of the recirculation region is also affected prominently at this frequency, suggesting that mass depletion of the recirculation region may be due to shedding observed during this oscillatory cycle. The separation bubble exhibits wavelike behaviour in the $St_{L}\sim 0.5$ band, consistent with Kelvin–Helmholtz shedding, which is identified downstream of the interaction region in wall pressure spectra. However, the mass budget of the reverse flow region does not indicate significant dynamics at this frequency, and rather it is observed that the bubble simply rides the waves. While it is more challenging to draw conclusions at this higher frequency regarding how the Kelvin–Helmholtz shedding is influenced by linear mechanisms due to a greater degree of contamination from highfrequency turbulence, the results are not inconsistent with the possibility that this shedding is also sustained by linear mechanisms evident throughout the interaction.
Acknowledgements
The authors acknowledge several insightful conversations with Dr P. Schmid, Dr L. Agostini and Dr S. Unnikrishnan, and appreciate the sponsorship of the US Air Force Office of Scientific Research (Monitor: Dr I. Leyva, contract FA95501410167). M.C.A. is supported through the US National Defense Science and Engineering Graduate Fellowship (NDSEG) Program. The simulations described were performed with grants of computer time from the US Department of Defense HPCMP shared resources at AFRL, NAVO and ERDC. Some figures have been made with complimentary licenses of FieldView obtained from Intelligent Light under the University Partners Program.