1. Introduction
Proposed next-generation supersonic aircraft consider progressively complex engine configurations to achieve demanding flight requirements, motivating the need to understand and control ever more intricate jet flow physics. A modern variable-cycle engine design, of increasing interest over the last decade, is discussed by Simmons (Reference Simmons2009) and shown schematically in figure 1. In addition to a rectangular cross-section to facilitate easier airframe integration and drag reduction, the arrangement utilizes a single-sided expansion ramp nozzle (SERN) with a unique independently modulated, cooled bypass stream to continuously match inlet airflow and engine demand while simultaneously serving as a cool heat sink to dissipate aircraft heat loads. Under suitable conditions, the addition of such streams has also been shown to reduce noise (Papamoschou & Debiasi Reference Papamoschou and Debiasi2001; Berry et al. Reference Berry, Magstadt, Glauser, Ruscher, Gogineni and Kiel2016). Numerous variations of the multistream nozzle design have been previously studied by Hromisin (Reference Hromisin2019), Papamoschou et al. (Reference Papamoschou, Phong, Xiong and Liu2016), Hromisin et al. (Reference Hromisin, Akatsuka, McLaughlin and Morris2019), Tinney, Panickar & Vogel (Reference Tinney, Panickar and Vogel2018) and Ruscher et al. (Reference Ruscher, Gogineni, Kiel, Giese and Viswanath2016) due to the advantages provided by the tertiary bypass stream and promising noise reductions, making this a highly coveted design for implementation into future supersonic air vehicles.
Illustration of generic three-stream engine (Simmons Reference Simmons2009).
The present study analyses the supersonic flow produced at the exit of the engine shown in figure 1, where it is assumed that the indicated core and fan streams are fully mixed before exiting the nozzle. A cutaway of the nozzle cross-section region of interest is shown in figure 2(a). Thus, the nozzle contains two canonical flows: a main core stream and a bypass stream with a splitter plate separating the two, while an aft-deck is used to represent the structure of an aircraft wing or fuselage. The cross-section becomes complete with the addition of the upper SERN nozzle surface.
(a) Nozzle configuration and (b) nomenclature.
The aerodynamic flow at the engine exit contains various turbulent mixing layers, which have been previously studied by Delville et al. (Reference Delville, Ukeiley, Cordier, Bonnet and Glauser1999) and Ukeiley et al. (Reference Ukeiley, Cordier, Manceau, Delville, Glauser and Bonnet2001). Supersonic jets have been known to produce jet noise that is generally associated with turbulent mixing, broadband shock-associated noise and screech tones. Investigation into jet noise was first pioneered by Lighthill (Reference Lighthill1952, Reference Lighthill1954) and later extended to supersonic turbulent shear layers by Phillips (Reference Phillips1960), who predicted the sound field radiates as eddy Mach waves that are produced along the shear zone. Moreover, interactions between the shock cells and shear layer instabilities downstream of the nozzle exit generate acoustic waves that are mostly broadband. In the scenario where the upstream propagating sound waves excite the shear layer near the nozzle lip, screech tones are produced from the closure of a feedback loop as Powell (Reference Powell1953) first suggested. Under certain conditions, harmonics of the screech tone can also be observed (Tam, Parrish & Viswanathan Reference Tam, Parrish and Viswanathan2014).
In the present work, the nozzle houses a supersonic shear layer that is generated from the two incoming streams. The flow downstream of the relatively thick splitter plate is characterized by undesirable flow features whose influence permeates other regions in the entire cross-section of the nozzle exit and also the plume, in the form of tonal noise and unsteady surface loading. The main flow features are vortex shedding (VS) and its consequences, including the formation of a relatively pronounced unsteady shock system that traverses the flow, interacts with and reflects from the SERN back on the deck surface. Extensive computational and experimental work has been performed in an attempt to control the flow. Because the dominant flow physics of shedding vortices from the splitter plate contributes most to the far-field resonant tone, these efforts have been directed towards perturbing the flow near the splitter plate trailing edge (SPTE). Passive control in the form of a spanwise wavy splitter plate was investigated experimentally by Gist (Reference Gist2022) and computationally by Doshi et al. (Reference Doshi, Gaitonde, Gist and Glauser2022). A key finding was that the streamwise vorticity induced by the wavy splitter plate broke down the coherent structures associated with the signature high-frequency tone found in far-field acoustics, resulting in a reduction of the dominant tone.
The success of this passive control technique prompted investigation into active control strategies such that similar control outcomes could be achieved, allowing for the potential development of an efficient and flexible sensor–actuator control system capable of operating under various flight conditions and with minimal energy input. Recently, Kelly & Glauser (Reference Kelly and Glauser2023) implemented microjet arrays that blow steady, sonic air into the flow. Because of experimental limitations, the actuators were placed on the aft-deck, beneath the SPTE, rather than on the potentially more effective locations, such as the splitter plate itself. Results with blowing at different angles indicated encouraging control authority on the shock system, which only modestly altered the tone. The present work lifts experimental constraints by performing a parametric study examining how actuator placement at different locations on the splitter plate affects the resonant tone and other flow characteristics.
Unsteady simulations using the Navier–Stokes equations are performed to obtain the data for the baseline and controlled conditions. An understanding of the simulated flow fields may be derived from modal analyses of the unsteady data. Various data-driven modal analysis techniques are appealing for their ability to decompose a turbulent flow field into modes that allow one to identify physically important flow features. Taira et al. (Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017), Taira et al. (Reference Taira, Hemati, Brunton, Sun, Duraisamy, Bagheri, Dawson and Yeh2020) and Rowley & Dawson (Reference Rowley and Dawson2017) have presented an overview and framework for several well-known modal decomposition methods and their applications. Proper orthogonal decomposition (POD) is a popular technique first introduced by Lumley (Reference Lumley1967) as a means to extract coherent structures from turbulent flow fields. It seeks an optimal set of basis functions that represent the given dataset, where each POD mode captures a portion of the total energy. In the frequency-based form of POD, known as the spectral proper orthogonal decomposition (SPOD) (Glauser, Leib & George Reference Glauser, Leib and George1987), dominant structures are instead extracted across a range of frequencies. It has proven to be a powerful tool in characterizing a turbulent flow and has seen increasing use across the fluid dynamics community. The properties of SPOD and its relationship to other approaches have been elaborated by Towne, Schmidt & Colonius (Reference Towne, Schmidt and Colonius2018). Examples, where SPOD variants have been leveraged to analyse both experimental and numerical datasets, include turbulent jets (Glauser et al. Reference Glauser, Leib and George1987; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Akamine et al. Reference Akamine, Tsutsumi, Okamoto, Teramoto and Nonaka2022), turbulent wakes (Araya, Colonius & Dabiri Reference Araya, Colonius and Dabiri2017; Nidhan et al. Reference Nidhan, Chongsiripinyo, Schmidt and Sarkar2020) and flows past airfoils (Abreu, Cavalieri & Wolf Reference Abreu, Cavalieri and Wolf2017; Ribeiro & Wolf Reference Ribeiro and Wolf2017). In the present study, the SPOD is applied to numerical datasets obtained from direct numerical simulations to identify flow structures at each frequency.
To inform actuator placement from a sensitivity and receptivity perspective and to characterize the effects of control, the simulated results are subjected to the resolvent analysis, which yields forcing and response dynamics of the simulated flow. The resolvent analysis is an operator-based modal technique to extract the input–output dynamics of a linear system. It may be viewed as complementary to global linear stability analysis which predicts the asymptotic evolution of disturbances in a base flow through an eigenvalue problem (McKeon & Sharma Reference McKeon and Sharma2010; Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017). In contrast, the resolvent analysis predicts forcing and response modes for chosen wavenumbers and frequencies through the singular value decomposition of the resolvent operator, which gives optimal energy amplifications between forcing and response modes. This yields physical insights into control designs, and is particularly beneficial when searching for effective active control configurations because the parameter space is extensive and costly to test with resolved simulations. Moreover, this approach facilitates investigation of the input–output dynamics of a system under the influence of control.
The paper is organized as follows. The active control configuration, computational approach and methods used to analyse the supersonic flow generated from a SERN are detailed in § 2. Results are presented in § 3, where instantaneous and mean flow modifications are discussed along with the overall instabilities of each case and findings from various modal analyses. Finally, concluding remarks are given in § 4.
2. Methods and set-up
2.1. Physical model problem
The supersonic nozzle flow at the exit of the engine of figure 1 is modelled as shown in figure 2, and includes the splitter plate with the aft-deck below and the SERN surface above. The nozzle dimensions are summarized in table 1. All length scales in the subsequent discussion are non-dimensionalized by the splitter plate width 
$W$, where 
$W$ corresponds to an experimental length scale of 82.296 mm. Further details on the geometry information may be found in the works by Simmons (Reference Simmons2009), Berry (Reference Berry2016), Magstadt (Reference Magstadt2017) and Stack (Reference Stack2019).
Nozzle dimensions (according to figure 2b) normalized by the splitter plate width 
$W$.
As noted earlier, among the three streams produced by the engine, the core and fan streams (figure 1) are well mixed before entering the SERN, forming one uniform main stream entering the nozzle (figure 2a). In the numerical set-up, only two physical streams are simulated for the nozzle flow: a supersonic main stream and a sonic bypass stream, which are separated by a splitter plate. The two streams mix before passing through and exiting the nozzle. The operating conditions of each stream are the same as those of prior experimental and numerical studies by Berry (Reference Berry2016), Magstadt (Reference Magstadt2017) and Stack (Reference Stack2019). The engine core flow operates at a nozzle pressure ratio (
${\rm NPR}= P_{total}/P_{ref}$) of 
${\rm NPR}_{core} = 4.25$, while the bypass stream operates at 
${\rm NPR}_{bypass} = 1.89$, corresponding to Mach numbers of 1.6 and 1.0, respectively, based on the isentropic relation. The nozzle temperature ratio (
${\rm NTR}= T_{total}/T_{ref}$) is set to unity for both flows, and is representative of the unheated jet utilized in the experiment by Berry (Reference Berry2016) and Magstadt (Reference Magstadt2017).
2.2. Active flow control configuration
Active flow control in the form of microjet actuation is introduced at different locations on the splitter plate near its trailing edge. Jets in crossflow have been extensively studied (Cortelezzi & Karagozian Reference Cortelezzi and Karagozian2001; Karagozian Reference Karagozian2014; Liu et al. Reference Liu, Sun, Liang, Cai and Wang2018; Li et al. Reference Li, Sun, Yu, Liang, Liu, Zhao and Huang2020) and used in active control efforts in a multitude of supersonic fluid flow problems (Karagozian Reference Karagozian2010; Guo, Liu & Zhang Reference Guo, Liu and Zhang2017; Guo et al. Reference Guo, Lin, Bu, Bai and Chao2018; Liu, Zhang & Liu Reference Liu, Zhang and Liu2019). The specific control objectives of the present work are to suppress the resonant tone, attenuate shock strength and reduce the surface loading on the aft-deck caused by the intense flow unsteadiness. Figure 3 depicts the details of the control configuration. Three locations are considered independently, in which the actuator is placed on either the splitter plate top surface (
$\text {SP}_{T}$), or the splitter plate bottom surface (
$\text {SP}_{B}$) or the vertical surface of the splitter plate trailing edge (
$\text {SP}_{TE}$). These choices are motivated by the resolvent analysis, discussed below, which indicated that this region is most receptive to external disturbances.
Schematic of the active flow control configuration. Each slotted actuation surface is indicated with a thick red line. Here 
$\text {SP}_{T}$, splitter plate top surface; 
$\text {SP}_{TE}$, SPTE surface; 
$\text {SP}_{B}$, splitter plate bottom surface.
The microjet actuator is characterized by two parameters: the location and angle of injection, 
$\psi$. At each location, steady blowing is introduced at different injection angles over a range of 
$\psi \in [-60^\circ, 90^\circ ]$, where 
$\psi$ is measured anticlockwise from the streamwise direction, 
$x$. The momentum coefficient, 
$C_\mu$, and the jet-to-mainstream momentum flux ratio, 
$J$, are defined as
\begin{equation} C_\mu = \frac{\dot m_{a} V_{a}}{\dfrac{1}{2}\bar{\rho}_{e} \bar{U}_{e}^2 A_{e}}\quad \text{and}\quad J = \frac{\rho_{a}V_{a}^2}{\bar{\rho}_{e}\bar{U}_{e}^2}, \end{equation}
where 
$\dot m_{a}$ is the actuation's blowing mass flow rate, 
$V_{a}$ is the centre local jet velocity, the sign of 
$\overline {({\cdot })}$ denotes spatial- and time-averaged quantity, 
$\rho _{e}$ is the density at the nozzle exit, 
$U_{e}$ is the velocity at the nozzle exit, 
$A_{e}$ is the area of the nozzle exit and 
$\rho _{a}$ is the density of the microjet. A hyperbolic-tangent velocity profile is prescribed for all the control cases, as shown in figure 3, with a centre local Mach number of 1.0 at the slot. The actuator location, injection angle, momentum coefficient and jet-to-mainstream momentum flux ratio of all the control configurations considered in the present work are tabulated in table 2.
Summary of the active control configurations: 
$\text {SP}_{T}$, splitter plate top; 
$\text {SP}_{TE}$, trailing edge; 
$\text {SP}_{B}$, bottom surfaces; 
$\psi$, injection angle; 
$\Delta x/W$, streamwise distance from actuation centre to SPTE; 
$C_\mu$, momentum coefficient; 
$J$, jet-to-mainstream momentum flux ratio.
2.3. Numerical model
In prior numerical studies, the dominant flow physics within the nozzle, except very near the sidewalls, were found to be largely two-dimensional (2-D). Figure 4 depicts the three-dimensional (3-D) simulations performed by Doshi (Reference Doshi2023) and Stack (Reference Stack2019), where figure 4(a,b) illustrate the shock topology of the full SERN configuration. Within the nozzle, the shock structures are observed to be largely 2-D, while shocks originating from the corner regions at the nozzle exit result in the formation of a 3-D shock train. Since the high-frequency tone associated with this jet was found to originate from the mixing of the two incoming flows within the nozzle, a 3-D simulation of the isolated shear layer was performed and is shown in figure 4(c,d). It was found that the coherent structures associated with the shedding instability were also primarily 2-D. Since we seek to perform a relatively large number of parametric studies, a 2-D rectangular domain is used to simulate the flow.
The 3-D simulation results by Doshi (Reference Doshi2023) and Stack (Reference Stack2019). Isosurfaces of (a) 
$\boldsymbol {U} \boldsymbol{\cdot} \boldsymbol {\nabla } P = 0.3$ and (b) 
$\boldsymbol {U} \boldsymbol{\cdot} \boldsymbol {\nabla } P = 0.2$ for the full SERN configuration. (c,d) Instantaneous visualization of the 
$Q$-criterion for an isolated shear layer.
Two-dimensional direct numerical simulations are performed using the solver, CharLES (Brès et al. Reference Brès, Jordan, Jaunet, Le Rallic, Cavalieri, Towne, Lele, Colonius and Schmidt2018; Brès & Lele Reference Brès and Lele2019). Due to this, the subsequent discussion is limited to the 2-D turbulent structures in the flow. The solver uses a second-order finite-volume method and a third-order Runge–Kutta temporal scheme to solve the compressible Navier–Stokes equations, and the relative-solution-ENO scheme is utilized to capture shocks. The unsteady compressible Navier–Stokes equations are non-dimensionalized using the reference values of 
$\rho _{ref} = 1.173$ kg m
$^{-3}$, 
$c_{ref} = 347.189$ m s
$^{-1}$, 
$P_{ref} = \rho _{ref} c_{ref}^2$ and 
$T_{ref} = 300$ K. The Reynolds number is of 
$Re = \rho _\infty U_{jet} D_h/\mu _\infty = 1.5\times 10^5$, where 
$\rho _\infty$ is the free-stream density, 
$U_{jet}$ is the nozzle exit velocity calculated based on the isentropic relation and 
${\rm NPR}_{core} = 4.25$, 
$D_h$ is the hydraulic diameter and 
$\mu _\infty$ is the free-stream dynamic viscosity. The hydraulic diameter, 
$D_h$, is calculated based on the rectangular nozzle exit, defined as
\begin{equation} D_h = \frac{2h_{NE}W}{h_{NE}+W}, \end{equation}
where 
$h_{NE}$ and 
$W$ are the nozzle exit height and width, respectively. Here, the nozzle width is equivalent to the splitter plate width.
Figure 5(a) shows the computational domain and flow configuration used in the present study. The length and width of the domain are 13.253
$W$ and 11
$W$, respectively. The streamwise and wall-normal directions are denoted by 
$x$ and 
$y$ with corresponding velocity components of 
$u$ and 
$v$, where a Cartesian coordinate system is used. As shown in figure 5(b), a structured mesh with non-uniform spacing is used, where grids along no-slip walls and around primary jet flows are refined to resolve boundary layers and small-scale structures. A grid resolution study is performed for three mesh sizes containing 
$0.7 \times 10^6, 1.0~\times 10^6$ and 
$2.0~\times 10^6$ points, where each grid further refines the region surrounding the primary jet flow. The time-averaged velocity profiles at the nozzle exit, pressure spectra and turbulence intensities (TI) of the three shear layers are compared between each grid in figure 6. The TI is calculated using the velocity components as
\begin{equation} \text{TI}(t) = \sqrt{\frac{\dfrac{1}{2}(u'(t)^2 + v'(t)^2)} {\bar{u}^2 + \bar{v}^2}}, \end{equation}
where 
$u'(t)$ and 
$v'(t)$ are the time-varying fluctuating components of 
$u$ and 
$v$, respectively. While all grids are able to capture the resonant tone, the difference in TI between the 
$1.0 \times 10^6$ and 
$2.0 \times 10^6$ grids is small. Moreover, the corresponding velocity profiles agree well near the walls. The 
$1.0 \times 10^6$ mesh is chosen in favour of the computational cost and is used for the remainder of this study (figure 5b).
(a) Schematic of computational domain set-up and flow configuration. Nozzle geometry is shown in grey, and sponge layers are indicated in blue. (b) Grid topology with instantaneous vorticity.
Time-averaged velocity profiles at the nozzle exit, pressure spectra and TI for three grid resolutions: 0.7, 1.0 and 
$2.0\times 10^6$. The case with 
$2.0 \times 10^6$ cells is shown dashed for clarity in the smaller window.
A coflow of Mach number 
$M$ = 0.01 is employed at the left-hand boundary to mimic the experimental set-up used by Berry (Reference Berry2016) and Magstadt (Reference Magstadt2017) at the Skytop Turbulence Lab at Syracuse University. Characteristic boundary conditions are applied at the main and bypass stream boundaries. The inflow conditions are determined using the designed nozzle pressure ratio and nozzle temperature ratio for both streams, and are prescribed to be uniform laminar flows. All the walls of the SERN, splitter plate and aft-deck are assigned to be adiabatic with a no-slip condition. An upwind scheme (Brès et al. Reference Brès, Ham, Nichols and Lele2017) is employed in a region of length 8
$W$ and width 6
$W$. Sponge layers are placed outside of the upwind domain along the top, bottom and outflow boundaries. In this region, a source term is used to damp outgoing waves and prevent reflection from boundaries into the domain (figure 5a).
A comparison between the 3-D simulation performed by Stack (Reference Stack2019) and the current 2-D results is shown in figure 7. Similar shock train components as captured by Stack (Reference Stack2019) are observed in the 2-D simulation, most notably Shocks 1, 2A, 2B and 4. The effects of choices made to enable parametric studies have a modest effect on the solution. Thus, Shock 1 in the 3-D simulation impinges along the SERN closer to the nozzle lip since Stack's Reynolds number is one order of magnitude higher (
$Re_{D_h} = 2.70\times 10^6$) than the current 2-D study. In turn, Shock 2A develops farther upstream in the present work, and the reflected shocks from 2A and 2B do not coalesce into Shock 3, which also consists of the 3-D corner shocks, as shown in figure 4(b). However, Shocks 2B and 4 are observed in the same locations for both the 2-D and 3-D configurations. The lowering of Reynolds number and use of 2-D simulations substantially reduces the required computational resources without compromising the primary mechanisms of interest behind the development of Shocks 1 and 2A. In fact, the characteristic high-frequency tone captured in the present work is the same as that reported in simulation and experiment, as will be discussed later.
Comparison between the time-averaged 
$\boldsymbol {U}\boldsymbol{\cdot} \boldsymbol {\nabla } P$ from the 3-D simulation by (a) Stack (Reference Stack2019) and (b) 2-D results in the present work.
2.4. Spectral proper orthogonal decomposition
To gain insight into the complex flow dynamics, the SPOD algorithm, a method initiated by Glauser et al. (Reference Glauser, Leib and George1987) and analysed recently by Towne et al. (Reference Towne, Schmidt and Colonius2018), is applied on the baseline and representative control cases to extract coherent flow structures associated with prominent frequencies. In the SPOD algorithm shown in figure 8, an instantaneous state of the flow 
$\boldsymbol {q}(\boldsymbol {x},t)$ at time 
$t_k$ can be rearranged into a column vector 
$\boldsymbol {q}_k \in \mathbb {R}^N$, where 
$N$ is equal to the number of data points times the number of flow variables. Given 
$M$ equally spaced snapshots, a large data matrix 
$\boldsymbol {Q}$ is constructed such that 
$\boldsymbol {Q} = [\boldsymbol {q}_1\ \boldsymbol {q}_2 \ \cdots \ \boldsymbol {q}_M] \in \mathbb {R}^{N \times M}$. This snapshot matrix is then partitioned into overlapping blocks 
$\boldsymbol {Q}^{(n)} = [\boldsymbol {q}_1^{(n)} \ \boldsymbol {q}_2^{(n)} \ \cdots\ \boldsymbol {q}_{N_{FFT}}^{(n)}] \in \mathbb {R}^{N \times N_{FFT}}$ with each block containing 
$N_{FFT}$ snapshots. A discrete Fourier transform is then applied on each block to obtain 
$\hat {\boldsymbol {Q}}^{(n)} = [\hat {\boldsymbol {q}}_1^{(n)} \ \hat {\boldsymbol {q}}_2^{(n)} \ \cdots\ \hat {\boldsymbol {q}}_{N_{{FFT}}}^{(n)}]$ with 
$\hat {\boldsymbol {q}}_k^{(n)}$ being the Fourier component at frequency 
$f_k$ in the 
$n$th block. The Fourier coefficients at each frequency 
$f_k$ from each block are rearranged into a new data matrix 
$\hat {\boldsymbol {Q}}_{f_k} = \sqrt {\kappa }[\hat {\boldsymbol {q}}_k^{(1)} \ \hat {\boldsymbol {q}}_k^{(2)} \ \cdots\ \hat {\boldsymbol {q}}_{k}^{(N_{{blk}})}]$ where 
$\kappa = \Delta t/(sN_{{blk}})$. Here, 
$N_{{blk}}$ is the number of blocks and 
$s = \sum _{j=1}^{N_{{FFT}}} w_j^2$, where the scalar weights 
$w_j$ can be used to reduce spectral leakage due to the non-periodicity in each block. The cross-spectral density tensor at frequency 
$f_k$ can then be estimated as 
$\boldsymbol {S}_{f_k} = \hat {\boldsymbol {Q}}_{f_k} \hat {\boldsymbol {Q}}_{f_k}^*$. An eigenvalue decomposition is then performed on an 
$N \times N$ matrix,
\begin{equation} \boldsymbol{S}_{f_k} \boldsymbol{\varGamma \varPhi}_{f_k} = \boldsymbol{\varPhi}_{f_k} \boldsymbol{\varLambda}_{f_k}, \end{equation}
with 
$\boldsymbol {\varPhi }_{f_k}$ being the SPOD modes that are ranked by the corresponding eigenvalues given by the diagonal matrix 
$\boldsymbol {\varLambda }_{f_k}$. It is noted that the choice of norm is expressed through the discrete weight matrix 
$\boldsymbol {\varGamma }$. A detailed guide to the SPOD algorithm can be found in the works by Towne et al. (Reference Towne, Schmidt and Colonius2018) and Schmidt & Colonius (Reference Schmidt and Colonius2020).
Schematic of SPOD algorithm. Each slice in 
$\boldsymbol {Q}$ represents a snapshot in time; DFT, discrete Fourier transform; CSD, cross-spectral density tensor.
In the present work, SPOD is computed using the density variable, where the simulation results are interpolated onto a 
$500\times 250$ uniform grid prior to performing SPOD. A grid study has been performed to ensure the SPOD results were not grid dependent. The weight matrix 
$\boldsymbol {\varGamma }$ is assigned to be an identity, which directly measures the variance of the data and is the natural choice of the norm if all data points are considered to be equally significant (Schmidt & Colonius Reference Schmidt and Colonius2020). For all the SPOD analysis performed in the present work, 5000 consecutive snapshots covering a time period of 
$t \overline {U_e}/D = 2.07$ are used, where 
$D$ is the streamwise length of the computational domain shown in figure 5. Each snapshot is separated by 
$\Delta t \overline {U_e}/D = 4.14 \times 10^{-4}$. The snapshot matrix 
$\boldsymbol {Q}$ is partitioned into 
$N_{{blk}} = 15$ blocks with 75
$\%$ overlap such that each block contains 
$N_{{FFT}} = 1111$ snapshots.
2.5. Resolvent formulation
To gain physical insights into active control configurations, resolvent analysis is performed using the turbulent mean flow to find optimal energy amplifications. Comprehensive details of the resolvent formulation can be found in the works by McKeon & Sharma (Reference McKeon and Sharma2010), Sun et al. (Reference Sun, Liu, Cattafesta, Ukeiley and Taira2020) and Yeh & Taira (Reference Yeh and Taira2019). Here, the compressible Navier–Stokes equations are expressed as
\begin{equation} \frac{\partial \boldsymbol{q}}{\partial t} = \mathcal{N}(\boldsymbol{q}), \end{equation}
where 
$\mathcal {N}$ is the nonlinear Navier–Stokes operator applied on the state variables 
$\boldsymbol {q} = [\rho,u,v,w,T]^{\rm T}$. Here 
$\rho$ is the density, 
$u, v, w$ are the velocity components, and 
$T$ is the temperature. The state variables can be Reynolds decomposed into 
$\boldsymbol {q}(x,y,z,t) = \bar {\boldsymbol {q}}(x,y) + \boldsymbol {q}'(x,y,z,t)$, where 
$\bar {\boldsymbol {q}}$ and 
$\boldsymbol {q}'$ are the mean and fluctuating components, respectively. Substituting the decomposed state variables into (2.5) yields
\begin{equation} \frac{\partial \boldsymbol{q}'}{\partial t} = \mathcal{L}_{\bar{\boldsymbol{q}}}(\boldsymbol{q}') + \boldsymbol{f}', \end{equation}
where 
$\mathcal {L}_{\bar {\boldsymbol {q}}}$ is the linearized Navier–Stokes operator about a base state and 
$\boldsymbol {f}$ denotes higher-order nonlinear terms or external forcing. In the present work, the operator 
$\mathcal {L}_{\bar {\boldsymbol {q}}}$ is constructed using a time-averaged flow. The terms 
$\boldsymbol {q}'$ and 
$\boldsymbol {f}'$ are then represented as a sum of Fourier modes with
\begin{equation} \boldsymbol{q}'(x,y,z,t) = \hat{\boldsymbol{q}}_{\omega,\beta}(x,y)\ {\rm e}^{{\rm i}(\beta z-\omega t)} \end{equation}and
\begin{equation} \boldsymbol{f}'(x,y,z,t) = \hat{\boldsymbol{f}}_{\omega,\beta}(x,y)\ {\rm e}^{{\rm i}(\beta z-\omega t)}, \end{equation}
where 
$\omega$ is a complex-valued frequency and 
$\beta$ is a real-valued spanwise wavenumber. Substituting the Fourier expansions into (2.6) brings the linearized Navier–Stokes equations to Fourier space as
\begin{equation} {-}{\rm i}\omega \hat{\boldsymbol{q}}_{\omega,\beta} = \mathcal{L}_{\bar{\boldsymbol{q}}}(\hat{\boldsymbol{q}}_{\omega,\beta} ;\beta) + \hat{\boldsymbol{f}}_{\omega,\beta}, \end{equation}where (2.9) can be rearranged and expressed as
\begin{equation} \hat{\boldsymbol{q}}_{\omega,\beta} = [-{\rm i}\omega \boldsymbol{I} - \mathcal{L}_{\bar{\boldsymbol{q}}}(\beta)]^{-1} ~ \hat{\boldsymbol{f}}_{\omega,\beta}. \end{equation}Here,
\begin{equation} \mathcal{H}_{\bar{\boldsymbol{q}}}(\omega,\beta) = [-{\rm i}\omega \boldsymbol{I} - \mathcal{L}_{\bar{\boldsymbol{q}}}(\beta)]^{-1} \end{equation}
is referred to as the resolvent operator and acts as a transfer function between the input 
$\hat {\boldsymbol {f}}$ and output 
$\hat {\boldsymbol {q}}$. Performing a singular value decomposition on the operator 
$\mathcal {H}_{\bar {\boldsymbol {q}}}(\omega,\beta )$ yields optimal forcing–response mode pairs with an associated gain, evaluated as
\begin{equation} \mathcal{H}_{\bar{\boldsymbol{q}}}(\omega,\beta) = \boldsymbol{\mathcal{Q}} \varSigma \boldsymbol{\mathcal{F}}^*, \end{equation}
where the singular vectors 
$\boldsymbol {\mathcal {Q}}$ and 
$\boldsymbol {\mathcal {F}}$ represent the response and forcing modes, respectively. The corresponding gains for each forcing–response pair are ranked in terms of their magnitude as 
$\varSigma = \text {diag}(\sigma _1,\sigma _2,\ldots,\sigma _m$). In the present study, the spanwise wavenumber is always chosen to be 
$\beta = 0$ such that focus is placed on the 2-D mechanisms, and Chu's norm for compressible flows (Chu Reference Chu1965) is used for the resolvent analysis (Towne et al. Reference Towne, Schmidt and Colonius2018; Yeh & Taira Reference Yeh and Taira2019).
In the event of an unstable base state, the stability of the linearized Navier–Stokes operator 
$\mathcal {L}_{\bar {\boldsymbol {q}}}$ is analysed as a precursor to the resolvent analysis (Sun et al. Reference Sun, Taira, Cattafesta and Ukeiley2017, Reference Sun, Liu, Cattafesta, Ukeiley and Taira2020; Liu et al. Reference Liu, Sun, Yeh, Ukeiley, Cattafesta and Taira2021). An eigenvalue decomposition of the linear operator is performed as
\begin{equation} \omega \hat{\boldsymbol{q}}_{\omega,\beta} = {\rm i}\mathcal{L}_{\bar{\boldsymbol{q}}}(\beta) \hat{\boldsymbol{q}}_{\omega,\beta}, \end{equation}
where 
$\omega$ is a complex-valued eigenvalue with the real 
$\omega _r$ and the imaginary 
$\omega _i$ components representing the frequency and growth rate, respectively. The most unstable eigenvalue of the system is utilized to determine a real-valued parameter 
$\alpha$, where 
$\alpha$ is chosen such that 
$\alpha > \max (\omega _i)$. This parameter introduces an exponential discount to the system, forming the discounted operator presented by Jovanovic (Reference Jovanovic2004). This limits the flow response to a finite time window characterized by 
$1/\alpha$, and focus is placed on the short-term dynamics. Introducing the discount parameter 
$\alpha$ into the original resolvent operator presented in (2.11) yields
\begin{align} \mathcal{H}_{\bar{\boldsymbol{q}},\alpha}(\omega,\beta) & = [-{\rm i}\omega \boldsymbol{I} - \{\mathcal{L}_{\bar{\boldsymbol{q}}}(\beta) - \alpha\boldsymbol{I}\}]^{-1}\nonumber\\ & = [-{\rm i}(\omega+{\rm i}\alpha)\boldsymbol{I} - \mathcal{L}_{\bar{\boldsymbol{q}}}(\beta)]^{-1} \end{align}
and 
$\mathcal {H}_{\bar {\boldsymbol {q}},\alpha }$ is referred to as the discounted resolvent operator. It is noted that choosing 
$\alpha = 0$ recovers the original formulation shown in (2.11), which represents a scenario in which an infinite-long time window is considered.
3. Results
3.1. Baseline flow features
In this section, the instantaneous and mean flow features of the baseline case are discussed. Spectral analysis is performed to identify and extract coherent structures that correspond to prominent frequencies. The dominant instabilities of the flow are investigated, and the discounted resolvent analysis is applied to gain insight into the forcing–response relations of the flow, such that optimal locations for introducing perturbations can be discovered.
3.1.1. Instantaneous and mean flow features
An instantaneous and mean flow field of the baseline case is shown in figure 9. Figure 9(a) displays shaded colour contours of the streamwise velocity that are visualized together with black-white contours of the pressure field, while figure 9(b) shows the time-averaged density gradient in the 
$x$-direction so that negative and positive values highlight expansions and shocks, respectively, in the supersonic flow. As shown in figure 9(b), an expansion fan forms as the main stream turns clockwise around the SPTE. Following this, a strong primary shock (S1) is observed, whose structure is comprised of an oblique shock and vortex-induced compression waves. The former is formed as the main stream aligns with the horizontal direction downstream of the SPTE, while the latter is associated with VS generated due to the mixing between the main and bypass streams. These affect the properties of the initial turning-associated oblique shock to form the S1 shock. As shown in figure 9(a), mixing of the main-bypass streams additionally causes a small separation region after the SPTE (SR1). The shedding vortices generated by the mixing of the two streams impinge on the aft-deck and convect downstream, forming a splitter plate shear layer (SPSL). The S1 shock travels across the core flow and impinges on the expansion ramp of the nozzle, inducing a shock-boundary-layer interaction (SBLI) and causing a shock-induced flow separation (SR2). The reflected R1 shock then impinges on the aft-deck, where it encounters the incoming vortices, and a vortex-shock–boundary-layer interaction (V-SBLI) occurs. A shock is observed emanating from the nozzle lip due to the static pressure at the nozzle exit being lower than the ambient, indicative of the overexpanded nature of the baseline configuration. This is consistent with prior experimental and numerical studies (Berry Reference Berry2016; Magstadt Reference Magstadt2017; Stack Reference Stack2019). When the flow exits the nozzle, the upper shear layer (USL) and lower shear layer (LSL) form with the ambient.
Instantaneous (a) and mean (b) flow fields for baseline case.
3.1.2. Spectral analysis
The frequency spectrum is first examined using the power spectral density (PSD) of the pressure time series at various point probes. The PSD is calculated using Welch's method (Welch Reference Welch1967) with a Hanning window and 75 % overlap. The data is sampled covering the entire data length of 
$t_D = t\overline {U_{e}}/D > 10$ for all data sets, and the sampling frequency is 
$f_s = D/\Delta t\overline {U_{e}} = 1.0526\times 10^4$. Figure 10 presents a non-dimensional PSD 
$\text {PSD}^* = 10\log _{10} (P_{xx} U_{jet}/q_\infty ^2 D_h)$ for two point probes, denoted as M1 and FF0, as well as their corresponding locations in the flow field, where 
$P_{xx}$ is the autospectral density and 
$q_\infty$ is the free-stream dynamic pressure. The M1 probe 
$(x,y)/W = (4.017, 0.940)$ is placed in the splitter plate wake and directly captures the shedding frequency of Strouhal number 
$St = fD_h/\overline {U_{e}} = 3.28$ as well as its harmonic. This is the same tone captured numerically by Stack (Reference Stack2019) (33 kHz) and experimentally by Berry (Reference Berry2016) (34 kHz) of approximately 
$St = 3.3$. This resonant tone has been found to dominate the flow field as well as the experimental far-field acoustics (Berry Reference Berry2016). Farther downstream, the FF0 probe 
$(x,y)/W = (4.700, 1.300)$ is placed near the nozzle lip and is shown to capture the same resonant tone. The FF0 probe will be used for the remainder of this paper to allow comparison between various control cases, as the M1 probe proves unreliable when actuation is introduced in the 
$\text {SP}_{TE}$ location (see figure 3) because it is dominated by the signal of the steady microjet.
Power spectral density of the pressure signal from two point probes in the baseline flow. Probe locations are shown in the lower left-hand corner.
To gain further insight into the base flow, SPOD is applied to the density variable to extract coherent structures from the flow field, where the identity weight matrix directly measures the variance of the data. Figure 11(a) shows the energy distribution for the baseline flow where 
$\varPhi _{i = 1}$ is the leading mode containing the most energy with increasing 
$i$ indicated with an arrow. Figure 11(b) shows the leading mode shapes corresponding to prominent peaks in the energy spectra. Peaks can be seen in the energy spectrum at the dominant frequency of 
$St = 3.28$ and its harmonics, similar to the M1 pressure spectrum shown in figure 10. Visualization of the leading mode at these frequencies reveals mode structures in the splitter plate shear layer, which further reinforces the dominance of the VS in the overall flow. Although these modes visualize the density variation, the shape of the modal structures indicates the shedding vortical structures induce pressure waves that leave the nozzle exit.
(a) The SPOD energy spectra on the density variable for the baseline flow. (b) Leading modes at representative frequencies.
3.1.3. Stability and resolvent analysis
Stability analysis is performed to identify the dominant instability for the baseline flow, and also to determine an appropriate discount parameter 
$\alpha$ for the discounted resolvent analysis. Figure 12(a) shows the eigenspectrum of the baseline flow with the spanwise wavenumber 
$\beta = 0$, where eigenvalues corresponding to shock-pattern modes, VS modes, and USL and LSL modes are denoted with different symbols, while all stable modes are in grey. Representative unstable modes in each category with the highest growth rates are shown in figure 12(b). It is observed that instabilities corresponding to shock patterns (
) occur at low frequencies, and the leading eigenmode is a stationary shock mode with 
$St = 0.00$. The VS instability (
) dominates the midfrequency to high-frequency range, and the leading VS eigenmode occurs near the resonant frequency of 
$St = 3.28$. There is additionally a USL and an LSL instability (
), however, the associated growth rate is relatively small. Based on these results of unstable modes, the discount parameter is chosen to be 
$\alpha D_h/\overline {U_e} = 3.0$. This value ensures that all unstable eigenvalues in the baseline and control cases (discussed later) will be enclosed. Generally, choosing a larger discount parameter will result in a decrease in the gain magnitude, however, the mode shapes are mostly unaffected. Prior studies using the discounted resolvent analysis on a laminar airfoil separation (Yeh & Taira Reference Yeh and Taira2019) and a turbulent cavity (Liu et al. Reference Liu, Sun, Yeh, Ukeiley, Cattafesta and Taira2021) have shown that introducing discounting at larger values removes local spikes in the distribution attributable to subdominant and spurious eigenmodes near the neutral axis, while preserving the overall profile of the gain distribution.
(a) Eigenspectrum for baseline flow. (b) Representative unstable modes for 
$\hat {u}$.
To uncover how the flow responds to external forcing in a control effort, resolvent analysis is performed with 
$\beta = 0$ to predict locations that are most receptive to external perturbations. The discounted resolvent operator is used due to the instabilities present in the system, as shown in figure 12. Forcing and response mode pairs are obtained at each frequency with a corresponding gain 
$\sigma _i$. Figure 13(a) shows a plot of the first 
$\sigma _1$ and second 
$\sigma _2$ singular values across a range of frequencies. Figure 13(b) shows the optimal forcing and response mode pairs at representative frequencies, where yellow–green contour lines represent forcing modes and blue–red contour lines represent response modes. In the baseline case, two local peaks are seen with forcing and response modes distributed in the USL, denoted as Mode (1), and splitter plate shear layer, denoted as Mode (2).
(a) Plot of the gain 
$\sigma _i$ across a range of frequencies. (b) Optimal forcing and response modes for 
$\hat {u}$. Yellow–green contour lines represent forcing modes, and blue–red contour lines represent response modes.
In Mode (1), the spatial structures of the response mode are observed to originate from inside the nozzle, near the location at which the primary shock impinges on the nozzle wall and induces flow separation. This suggests that the USL response may be associated with the separation in SR2 (figure 9a), and reducing the size of the separation can influence the USL instability. On the other hand, the splitter-plate shear layer forcing–response modes, denoted as Mode (2), appear at a frequency of 
$St = 5.71$ rather than the resonant tone. This is likely due to the non-normal relation between the two most dominant VS eigenvalues (see figure 12) when performing the discounted resolvent analysis. In a non-normal system matrix, transient growth may be observed even if all eigenvalues lie in the stable plane, i.e. energy amplification is possible prior to an exponential decay as time approaches infinity (Schmid & Brandt Reference Schmid and Brandt2014). Since the eigenvalues inherently describe the time-asymptotic behaviour for non-normal systems, they therefore fail to capture processes that occur over finite time scales. In the context of the discounted resolvent, which focuses on short-term dynamics, the identified optimal response can thus occur at a forcing frequency that is far from an eigenvalue, known as pseudoresonance (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993). Schmid & Brandt (Reference Schmid and Brandt2014) have suggested that the difference between the resolvent norm and the eigenvalue-based response must be attributed to the non-orthogonality of the eigenvectors.
Since the shedding instability is most prevalent in the flow and is the source of the primary resonant tone, the spatial location highlighted by Mode (2) is utilized to inform the actuator placement rather than Mode (1). The Mode (2) forcing is observed near the SPTE, indicating that introducing perturbations in this region has the most potential for influencing the development of the shear layer instability. Moreover, the leading gain 
$\sigma _1$ associated with Mode (2) is greater than Mode (1), suggesting that acting on Mode (2) may have a greater capacity of altering the mean flow. This motivates the choice of spatial locations for active flow control previously presented in figure 3, where actuators are placed in various locations around this region.
It is also observed that when considering the suboptimal distribution 
$\sigma _2$, forcing–response mode pairs pertaining to a single flow feature will follow a continuous curve throughout the frequency range. At 
$St = 4.63$, there is an intersection between two continuous curves, where the USL response (Mode (1)) moves into the suboptimal distribution while the VS response (Mode (2)) moves into the optimal distribution.
3.2. Active control: splitter plate top surface actuation
In this section, modifications introduced from the splitter plate top surface steady-blowing microjet actuation are discussed, and new microjet-induced flow features are analysed. The SPOD and resolvent analysis are applied to each control case to investigate changes in the forcing–response dynamics of the flow.
3.2.1. Instantaneous and mean flow features
When steady blowing is introduced along the top surface of the splitter plate (
$\text {SP}_{T}$ in figure 3) at high actuation angles 
$\psi \in [60^\circ, 90^\circ ]$, the primary shock becomes visibly absent, and two unsteady sets of pressure waves emerge from the microjet, as shown in figure 14. The instantaneous density gradient shows that actuating at 
$\psi = 90^\circ$ causes the unsteadiness to propagate farther upstream into the nozzle, and as discussed later, the aft-deck surface loading is also increased the most. Due to this, the flow structures associated with 
$\psi = 60^\circ$ actuation are chosen for further analysis.
Instantaneous and mean flow fields for control cases with microactuation introduced at the 
$\text {SP}_{T}$ location with (a) 
$\psi = 60^\circ$ and (b) 
$\psi = 90^\circ$.
Figure 15(a) shows a closer view of the microjet at 
$\psi = 60^\circ$ with isolines of Mach 1.0 and figure 15(b) illustrates the formation mechanism of the two sets of pressure waves. The introduction of the microjet causes separation upstream of the actuator; the associated shear layer introduces a Kelvin–Helmholtz instability that generates the first set of pressure waves as the instability grows. As the core flow is still expanding above the splitter plate and has yet to reach Mach 1.0, these waves are allowed to propagate upstream. When the unsteady instability reaches and interacts with the microjet, a second set of pressure waves, inclined at a higher angle relative to the first, is formed due to the flow-induced flapping motion of the actuation. The two sets of waves propagate upwards, and impinge and reflect off the SERN wall. This interaction between the main stream and the microjet causes the incoming main flow to be slightly unsteady as it mixes with the bypass stream, which inhibits the stacking of the compression waves induced by the shedding vortices after the splitter plate. Hence, the S1 shock no longer forms, which additionally eliminates the shock-induced separation region SR2 (figure 9).
(a) Close-up view with isolines of Mach 1.0 and (b) schematic of microjet proximity region for 
$\text {SP}_{T}$ at 
$\psi = 60^\circ$.
At low actuation angles 
$\psi \in [30^\circ, 45^\circ ]$, shown in figure 16, the separation upstream of the microjet is reduced, so that the Kelvin–Helmholtz instability is no longer present. When the main stream encounters the microjet, a new relatively steady and stronger shock is formed instead. At 
$\psi = 45^\circ$, this shock is weaker and inclined at a higher angle relative to the shock formed at 
$30^\circ$. In both these low actuation angle cases, the primary shock remains mostly unaffected and thus has minimal control authority over the shock train development.
Instantaneous and mean flow fields for control cases with microactuation introduced at the 
$\text {SP}_{T}$ location with (a) 
$\psi = 45^\circ$ and (b) 
$\psi = 30^\circ$.
3.2.2. Spectral analysis
The pressure spectra through the M1 and FF0 probes are shown for all actuation angles in figure 17. It is observed through the M1 location that the shedding frequency decreases at higher actuation angles, and at 
$60^\circ$, the amplitude of the resonant tone is significantly reduced. Meanwhile farther downstream, the FF0 location demonstrates the overall broadband content decreases with the angle. As high-angled actuation was found to be more beneficial towards the mitigation of the primary shock development, the most optimal angle at the 
$\text {SP}_{T}$ location is chosen to be 
$\psi = 60^\circ$.
The PSD of the pressure time series through the M1 and FF0 probes for all 
$\text {SP}_{T}$ control cases.
To further understand how the flow physics are modified, SPOD is performed on the density variable for the optimal case of 
$\psi = 60^\circ$. Figure 18(a) shows the SPOD energy spectra with the leading baseline distribution illustrated in black for comparison, and corresponding leading modes 
$\varPhi _1$ at representative frequencies are displayed in figure 18(b). The resonant tone corresponding to the VS instability is shifted to a lower frequency of 
$St = 2.48$ with a reduced amplitude, as was observed in the M1 pressure spectrum. This suggests a reduction in the energy content associated with the shedding instability when control is introduced in this configuration. An additional high-frequency peak at 
$St = 5.67$ is observed, where this frequency is associated with pressure waves originating from the microjet that propagate and reflect off the top surface of the SERN. Compared with the baseline flow, where the VS held the majority of the energy content, the 
$\text {SP}_{T} 60^\circ$ case disperses this energy into the unsteadiness generated from the actuator.
(a) The SPOD energy spectra on the density variable for the 
$\text {SP}_{T}$ optimal control case at 
$\psi = 60^\circ$. Leading baseline energy distribution illustrated for comparison. (b) Leading modes 
$\varPhi _1$ at the frequencies 
$St = 2.48$ and 
$St = 5.67$.
Although two separate sets of pressure waves can be distinguished in the instantaneous flow field, as shown earlier in figure 14(a), only one peak in the SPOD spectrum is found to correspond to the unsteady waves. To identify the frequencies corresponding to the two sets of pressure waves, the PSD of a probe, marked P1, and the SPOD of the streamwise density gradient is computed in the microjet proximity region. Figure 19(a) shows the location of the P1 probe 
$(x,y)/W = (3.966, 0.980)$ as well as the reduced spatial window used for SPOD, where P1 is placed directly behind the microjet within the unsteadiness to capture their signals. The results are then reinterpolated onto a 
$150 \times 100$ uniform grid within the reduced SPOD window to resolve the two sets of pressure waves. The SPOD of the instantaneous density field did not provide a clear delineation of the two sets of pressure waves of interest; however, the gradient of the density field successfully highlighted the mechanism of interest, as shown in figure 14, and is used instead.
Spectral analysis for the 
$\text {SP}_{T}$ control case at 
$\psi = 60^\circ$ near the microjet proximity region. (a) Location of P1 probe and reduced spatial window for SPOD. (b) The PSD of the pressure time series through P1. (c) The SPOD energy spectra on the streamwise density gradient variable. (d) Leading modes 
$\varPhi _1$ at representative frequencies.
Figure 19(b) shows the pressure spectrum through P1, where two frequencies 
$St = 1.24$ and 
$St = 4.88$ are captured. The SPOD is used as a tool to visualize the spatial structures at those frequencies, shown in figure 19(c,d). In the SPOD energy spectra, the same frequencies as those captured in the pressure spectrum are found, where the corresponding leading modes distinctly illustrate the two unsteady waves. The spatial structures of their modes indicate that the higher frequency 
$St = 4.87$ corresponds to the first set of waves, while the lower frequency 
$St = 1.24$ corresponds to the second set of waves. It is also noted that despite using a reduced spatial window, the frequencies 
$St = 2.48$ and 
$St = 5.67$ captured in the larger domain SPOD remain dominant in these energy spectra. Visualizing their modes additionally shows similar structures presented earlier in figure 18. Although the frequencies 
$St = 1.24$ and 
$St = 4.87$ directly correspond to the unsteady pressure waves, their combined influence appears through the 
$St = 5.67$ tone.
3.2.3. Instabilities and resolvent modes
Stability analysis is performed for all angles 
$\psi$ considered at the 
$\text {SP}_{T}$ location, where the linearized operator 
$\mathcal {L}_{\bar {\boldsymbol {q}}}$ is reconstructed using the corresponding time-averaged flow obtained from each control case. Figure 20 shows the eigenspectra for all 
$\text {SP}_{T}$ control cases, where only the unstable eigenvalues (
$\omega _i > 0$) are shown. At high angles 
$\psi \in [60^\circ, 90^\circ ]$, a new actuation-induced instability with a relatively high growth rate is found, denoted with (
). In this mode, the unsteadiness generated from the actuator can be observed, and the spatial structures of the VS instability develop from the microjet rather than the SPTE, as previously shown in figure 12(b). At 
$\psi = 90^\circ$, the new unsteady mode is found predominantly in the midfrequency to high-frequency range (
$St \in [3.30, 5.43])$, while decreasing the angle to 
$\psi = 60^\circ$ shifts this mode into a wider range of frequencies (
$St \in [0.81, 6.16])$. Although these eigenmodes occur with high growth rates compared with the baseline spectrum, the unsteadiness associated with microjet-induced pressure waves proves beneficial in eliminating the primary shock, as shown in figure 14. It is also noted that the new actuation-induced mode is more unstable at 
$\psi = 60^\circ$ compared with 
$\psi = 90^\circ$, suggesting this unstable behaviour prompted the 
$\psi = 60^\circ$ case to perform better than the 
$\psi = 90^\circ$ case.
(a) Eigenspectra for all control cases at the 
$\text {SP}_{T}$ location with a zoomed-in view shown in a dashed grey box. Stable modes are not plotted for clarity. (b) New actuation-induced unstable mode.
Additional shock and USL and LSL instabilities are also observed at high angles in figure 20(a); however, both their frequencies and growth rates are relatively small. As the actuation angle decreases to 
$\psi \in [30^\circ, 45^\circ ]$, the eigenspectrum becomes largely concentrated with shock and shear layer modes, where these instabilities are mostly found in the low-frequency range for 
$\psi = 45^\circ$ and in the midfrequency range for 
$\psi = 30^\circ$. The most unstable eigenmode across all control cases occurs at 
$\psi = 30^\circ$ and corresponds to a stationary shock instability, indicating that introducing actuation closer along the incoming flow direction amplifies the shock unsteadiness behaviour.
Discounted resolvent analysis is performed with 
$\alpha D_h/\overline {U_e} = 3.0$ for all control cases to assess how the forcing–response dynamics of the flow have been changed after actuation is introduced. Figure 21(a) shows the distribution of the leading gain 
$\sigma _1$ for all angles 
$\psi$, and is presented as a contour plot across a range of frequencies. In the baseline distribution, the location of the USL and VS forcing–response mode pairs, previously shown in figure 13, are denoted as (
) and (
), respectively. For each case, forcing and response modes are obtained and visualized at local frequency peaks and throughout the distribution. Figure 21(b) presents representative forcing (yellow–green contour lines) and response (blue–red contour lines) modes, where 
$\hat {u}$ is plotted.
(a) Distribution of the optimal resolvent gain 
$\sigma _1$ across a range of frequencies for all 
$\text {SP}_{T}$ cases. Symbols represent different types of forcing–response mode pairs. (b) Representative forcing and response modes with optimal gain for 
$\hat {u}$. Forcing modes coloured as yellow–green and response modes coloured as blue–red.
When control is introduced at 
$\psi = 90^\circ$, the VS mode is shifted to a lower frequency with a higher gain while the USL mode is shifted to a higher frequency with a low gain. At higher frequencies, a new forcing–response pair (
) appears with a lower energy amplification, with the spatial structure of this response mode evident in both the upper and splitter plate shear layers. The USL response appears at the highest frequency when 
$\psi = 60^\circ$, and the VS response can additionally be found at multiple frequencies. As 
$\psi$ decreases, the USL response moves towards a lower frequency while the VS response shifts to a higher frequency; at 
$\psi = 30^\circ$, the two types of forcing–response modes appear most similar to the baseline case. Introducing control at the splitter plate top surface generally increases the responsive frequency of the USL, while simultaneously lowering the frequency of the VS response. Further, the USL resolvent modes always appear as a broadband peak in the gain distribution, while VS modes appear as a sharp peak and can be distributed into multiple frequencies. Introducing control shifts the frequency at which these modes appear; with the specific values being sensitive to the angle of actuation.
Contrary to Mode (1) shown previously in figure 13, the spatial structures of all USL response modes in the control cases originate closer to the nozzle lip rather than from inside the nozzle, and the structures of the corresponding forcing modes are smaller. This suggests that the USL resolvent modes are spatially associated with the size of flow separation in SR2, as all 
$\text {SP}_{T}$ control cases either reduce or eliminate the size of separation (see figures 14 and 16). Similarly, the spatial structures of all VS response modes in the controlled flows are initiated from the top of the splitter plate surface, indicating this control location is influencing the spatial formation of the signature shedding instability such that it develops closer to the microjet. Additionally, as the actuation angle 
$\psi$ increases, the corresponding forcing modes become purely concentrated along the top surface of the splitter plate upstream of the actuator, rather than the SPTE proximity region as in the baseline case.
3.3. Active control: SPTE surface actuation
This section discusses the effects of actuation at the SPTE surface. In contrast to the splitter plate top surface actuation, this control location directly influences the size of the shedding vortices. The SPOD is applied to the most optimal case, while resolvent analysis is performed on all control cases.
3.3.1. Instantaneous and mean flow features
When steady blowing is introduced at the trailing edge surface of the splitter plate (
$\text {SP}_{TE}$ in figure 3), the microjet directly interferes with the mixing of the main and bypass streams and considerably smaller vortices are generated for all angles 
$\psi$ considered, as shown in figures 22(a) and 23. When the microjet is introduced in the same direction as the incoming main and bypass streams (
$\psi = 0^\circ )$, shown in figure 22(a), a secondary shock S2 is formed. Figure 22(b) shows a close-up view of the microjet at 
$\psi = 0^\circ$, and figure 22(c) illustrates the formation mechanism of the secondary shock. The microjet breaks down the original thick splitter plate shear layer into two thinner shear layers. As a result, upper and lower sets of vortical structures are generated, whose sizes are smaller compared with the vortices in the baseline flow (figure 9). Moreover, the size of the upper and lower sets of vortices relative to each other are disparate due to the difference in the local speed of the main and bypass streams when mixing with the microjet. Farther downstream of the SPTE, the upper and lower vortices interact and merge, essentially where the potential core due to the microjet collapses, forming a secondary shock S2.
(a) Instantaneous and mean flow fields for control cases at the 
$\text {SP}_{TE}$ location with 
$\psi = 0^\circ$. (b) Closer view and (c) schematic of microjet proximity region (not to scale).
Instantaneous and mean flow fields for control cases at the 
$\text {SP}_{TE}$ location with (a–c) 
$\psi > 0^\circ$ and (d–f) 
$\psi < 0^\circ$.
When 
$\psi$ is increased in the 
$30^\circ$ to 
$60^\circ$ range, shown in figure 23(a–c), the upper vortices become smaller, reducing the strength of the resulting compression waves and ultimately weakening the S1 shock. Eventually, S1 is formed almost entirely due to the flow deflection caused by the microjet, and the contribution from the upper vortex-induced compression waves is negligible. Meanwhile, the S1 shock is inclined at a higher angle relative to the baseline flow due to the microjet being angled upwards, and the S1 shock moves farther upstream as 
$\psi$ increases. Additionally, the size of the lower vortices relative to the upper vortices becomes larger, resulting in the strengthening of the S2 shock when the two sets of vortices interact and merge. This causes the separation in SR2 (figure 9) to be induced by the S2 shock rather than the S1 shock. At the optimal angle where 
$\psi = 30^\circ$, neither the S1 nor S2 shocks will induce the flow separation in SR2 near the nozzle lip. Simultaneously, the weakened S1 shock leads to a smaller difference in the streamwise velocity before and after the S1 shock, such that the flow processed by S2 is of higher speed. This results in the impingement of S2 farther downstream on the SERN and closer to the nozzle lip. As 
$\psi$ increases, the S1 and S2 shocks gradually separate and are formed apart from each other.
At negative angles 
$\psi < 0^\circ$, shown in figure 23(d–f), the upper vortices contribute weaker compression waves to the S1 shock. However, a stronger oblique shock is formed due to the microjet deflecting the main flow downwards before aligning with the horizontal direction. This leads to the overall strengthening of the S1 shock where further decreasing 
$\psi$ will increase the strength of S1 and consequently, the size of the separation in SR2.
3.3.2. Spectral analysis
The SPOD energy spectra and modes at representative frequencies are displayed in figure 24 for the optimal control case of 
$\psi = 30^\circ$. Peaks are captured in the leading energy distribution at multiple frequencies with reduced amplitudes compared with the baseline distribution. In the low-frequency range, this control strategy excites frequencies of 
$St = 1.33$ and 
$St = 2.92$, corresponding to the USL instabilities. Peaks are observed in the high-frequency range at 
$St = 7.18$ and 
$St = 21.80$. Modal structures at the former value show the lower shedding vortices generated from the microjet-bypass streams mixing, while the latter value highlights the upper shedding vortices generated from the mixing of the microjet and main streams. The higher gradients in the spatial structures of the upper shedding vortices additionally suggest they dissipate faster than the lower vortices.
(a) The SPOD energy spectra on the density variable for the 
$\text {SP}_{TE}$ optimal control case at 
$\psi = 30^\circ$. Leading baseline energy distribution illustrated for comparison. (b) Leading modes at representative frequencies.
3.3.3. Instabilities and resolvent modes
Eigenspectra obtained from the stability analysis of all control cases at the 
$\text {SP}_{TE}$ location are shown in figure 25. Similar to low angle actuation 
$\psi \in [30^\circ, 45^\circ ]$ at the 
$\text {SP}_{T}$ location, the shock (
) and USL and LSL instabilities (
) become significantly more unstable compared with the baseline spectrum. This might indicate that introducing actuation with a larger streamwise component will prompt an amplification in the shock and shear layer behaviour through mechanisms that agglomerate the flow deflection and vortex-induced compression waves that synthesize the primary S1 shock. These eigenmodes are distributed across a wider range of frequencies compared with the 
$\text {SP}_{T}$ low-angle actuation cases, shown in figure 20. The zoomed-in view of the dashed box region in figure 25 shows that for all angles considered, the USL instability (
) spans a large range of frequencies (
$St \in [0.0, 2.5]$) while the USL and LSL instabilities (
) appear together in a narrower and lower frequency range (
$St \in [0.0, 1.2]$). However, the leading eigenvalue with the lowest growth rate 
$\omega _i$ across the baseline and all control cases occurs at the optimal angle 
$\psi = 30^\circ$. This might suggest that reducing the size of the shedding vortices and achieving a balance between the primary and secondary shocks can most effectively suppress the baseline instabilities and also stabilize the controlled flow.
Eigenspectra for all control cases at the 
$\text {SP}_{TE}$ location with a zoomed-in view shown in a dashed grey box. Stable modes are not plotted for clarity.
Resolvent results of all control cases at the 
$\text {SP}_{TE}$ location are shown in figure 26 in the form of the leading gain 
$\sigma _1$ for all angles 
$\psi$ and representative forcing–response mode pairs. When actuation is introduced at 
$\psi \leq 0^\circ$, the USL and VS responses are both shifted to a higher frequency. At 
$\psi = -30^\circ$, the VS response is not present in the most optimal gain distribution 
$\sigma _1$, but rather in the suboptimal distribution 
$\sigma _2$. Meanwhile, a ‘hotspot’ corresponding to the USL response is evident around 
$St = 4.5$. This suggests that when actuation is introduced at a negative angle, the USL becomes progressively more responsive to perturbations, as indicated by the larger gain values. Simultaneously, the frequency of the USL response has minimal change across 
$\psi$, suggesting that although the USL becomes more responsive, it is not sensitive to changes in the actuation angle when 
$\psi \leq 0^\circ$.
(a) Contours of the optimal resolvent gain 
$\sigma _1$ across a range of frequencies for all 
$\text {SP}_{TE}$ control cases. (b) Representative forcing–response mode pairs with optimal gain for 
$\hat {v}$. Forcing modes coloured as yellow–green and response modes coloured as blue–red.
When 
$\psi \geq 0^\circ$, the corresponding gain and frequency of the USL response begins to decrease dramatically, implying this region of the flow becomes more stable; however, the USL becomes significantly more sensitive to the changes in 
$\psi$. The frequency of the VS response decreases as 
$\psi$ increases and begins to appear at multiple local peaks. This is likely due to the upper and lower shedding vortices (see figure 22c) becoming more apparent in the flow at positive actuation angles, whereas when 
$\psi < 0^\circ$, the lower vortices are inhibited in their development through the core flow (figure 23d–f) and the upper shedding vortices are largely dominant. Furthermore, the VS forcing modes concentrate closer to the SPTE, and the corresponding response modes propagate into the primary shock.
Compared with the baseline case, the USL response in all control cases originates closer to the nozzle lip, similar to the 
$\text {SP}_{T}$ control shown in figure 21. The size of the separation region SR2 is either reduced or eliminated in the controlled flows. This further signifies that the spatial response of the USL is largely associated with the SR2 separation. It is also observed that the USL resolvent modes consistently appear as a broadband peak while the VS modes appear as a sharp peak, as was observed in the 
$\text {SP}_{T}$ control.
3.4. Active control: splitter plate bottom surface actuation
When steady blowing is introduced at the bottom surface of the splitter plate (
$\text {SP}_{B}$ in figure 3), the bypass stream is obstructed for all angles considered and is most apparent in figure 27(a,b) when 
$\psi = -60^\circ$. The microjet stream reflects from the aft-deck, and the clockwise turn angle of the main stream around the SPTE is larger compared with the baseline flow. The main stream mixes with the microjet rather than the bypass stream, causing the flapping motion at the SPTE to be amplified. This generates stronger compression waves and strengthens the S1 shock. The interaction between the actuator and main stream additionally causes both the formation of the shedding vortices and the S1 shock to be spatially delayed such that both develop farther downstream from the SPTE compared with the baseline flow (figure 9). As a result, the S1 shock impinges on the SERN closer to the nozzle lip and induces a larger separation in SR2. Meanwhile, the deflected microjet from the aft-deck diverts the trajectory of the shedding vortices upward, causing the vortices to develop through the core flow where they follow a parabolic-like trajectory indicated by a black dashed line in figure 27(c).
Instantaneous flow fields for control cases at the 
$\text {SP}_{B}$ location: (a) closer view of microjet region; (b) 
$\psi = -60^\circ$; (c) 
$\psi = -45^\circ$; (d) 
$\psi = -30^\circ$.
The PSD at the FF0 probe near the nozzle lip (figure 28) shows that the resonant tone and all harmonics are shifted to a higher frequency. These frequency peaks are significantly amplified for all angles considered at this location. Generally, introducing steady-blowing control at the 
$\text {SP}_{B}$ location is ineffective in achieving the control objectives and will instead enhance the dominant tone and strengthen the primary shock.
Pressure spectra of the signal from the FF0 point probe for 
$\text {SP}_{B}$ actuation at the angles 
$-30^\circ$, 
$-45^\circ$ and 
$-60^\circ$.
3.5. Aft-deck surface loading
To assess the severity of the surface loading along the aft-deck, the integrated pressure fluctuation along the surface from the SPTE to the end of the aft-deck, reported as the root-mean-square after subtracting the mean, is plotted in figure 29(a) for all cases. Figure 29(b) shows the instantaneous density field of representative control cases to visualize changes in the trajectory and size of the shedding vortices. When steady blowing is introduced at the 
$\text {SP}_{T}$ location, the resulting surface loading along the aft-deck for all cases is within 
$\pm$9 % of the baseline state. At high angles where 
$\psi = 60^\circ$ and 
$90^\circ$, the upward flow deflection caused by the microjet produces slightly larger vortices, which increases the loading along the aft-deck. At low angles where 
$\psi = 30^\circ$ and 
$45^\circ$, the size of the vortices remains comparable to the baseline flow, resulting in minimal change in the surface loading.
(a) Integrated root-mean-square of pressure fluctuations along aft-deck for all control cases. (b) Instantaneous density field for representative control cases. Green arrows indicate microjet.
When actuation is introduced at the 
$\text {SP}_{TE}$ location, the vortex-induced surface loading along the aft-deck is reduced across all angles considered due to the weakened contact between the vortices and the aft-deck. The greatest reduction in surface loading occurs when 
$\psi < 0^\circ$. In the uncontrolled flow, the mixing of the main and bypass streams near the SPTE region results in the highest local pressure fluctuation on the aft-deck surface close to the SPTE. When microjet actuation is introduced at angles 
$\psi < 0^\circ$, the injected flow continuously impinges on the aft-deck near the SPTE region, leading to smaller local pressure fluctuations. Coupled with the smaller scale vortices, this yields the greatest reduction in the surface loading. As 
$\psi$ increases from 
$-30^\circ$ to 
$30^\circ$, the lower vortices (see figure 22c) that come into contact with the aft-deck become larger in size, leading to the increase in the surface loading. As 
$\psi$ increases to 
$60^\circ$, the surface loading begins to decrease. This is likely due to the lower vortices developing slightly above the aft-deck with reduced contact.
When actuation is introduced at the 
$\text {SP}_{B}$ location, the vortex trajectory is deflected upward, and the vortices convect downstream through the core flow. Since the vortices no longer develop along the aft-deck, there is a moderate reduction in the surface loading. However, this reduction is not as significant as introducing actuation at the 
$\text {SP}_{TE}$ location. This is because the size of the vortices remains comparable to the baseline flow, indicating that influencing the vortex size is most pivotal in reducing the aft-deck surface loading.
4. Conclusion
Control of a supersonic rectangular dual-stream jet flow with steady-blowing microjets is examined by parametrically varying spatial locations and injection angles in the vicinity of a splitter plate where the shedding instability due to mixing of the main and bypass streams produces a prominent undesirable resonant tone. In the baseline flow, the shedding vortices generate compression waves that coalesce to form an oblique shock to deflect the main flow. This relatively strong primary shock traverses through the core flow and impinges on the top of the nozzle wall, where the shock–boundary-layer interaction induces flow separation and affects the USL developing in the plume downstream. The shedding vortical structures additionally cause high surface loading along the aft-deck.
A combination of Navier–Stokes, resolvent, and stability analysis is used to explore mitigation and alleviation of the undesired baseline flow features. Discounted resolvent analysis on the turbulent mean baseline flow identifies optimal energy amplification and provides physical insight into effective active flow control configurations. The optimal forcing and response modes indicate that the region near the SPTE is most receptive to external perturbations, motivating the introduction of active control at the splitter plate top, trailing edge and bottom surfaces, respectively. The SPOD of representative control cases identify coherent structures in the flow field, while discounted resolvent analysis, repeated on all control cases, helps assess forcing–response dynamics changes due to actuation.
When the actuator is placed on the splitter plate top surface, the optimal injection angle 
$\psi$, measured from the horizontal, to attenuate the shock structures and dominant tone is observed to be 
$\psi = 60^\circ$. The success of this angle is traced to the introduction of a new instability upstream of the actuator, detected with SPOD as well as instability analysis of the controlled flow, which predicts a new eigenmode of high growth rate. The consequent microjet-induced waves prevent the main shedding vortex-induced compression waves from coalescing, and the primary shock then no longer forms. This has consequences on the far nozzle wall, where shock-induced flow separation is muted, and, for all injection angles, the resolvent response modes highlighting the USL originate closer to the nozzle lip. This indicates that the USL response is spatially associated with the shock-induced flow separation.
When actuation is introduced on the SPTE surface, the microjet breaks the original thick splitter plate shear layer into two thinner shear layers, giving rise to upper and lower sets of shedding vortices. The former generates weaker compression waves, reducing the primary shock strength. Farther downstream from the SPTE, the upper and lower sets of vortices interact and merge, forming a secondary shock. Depending on the actuation angle 
$\psi$, the flow separation near the nozzle lip is induced by either the primary or secondary shock; however, the 
$\psi = 30^\circ$ case achieves a balance between the strength of both shocks, such that neither induces separation. The SPOD analysis also captures the two sets of shedding vortices, which occur at high frequencies relative to the resonant tone found in the baseline flow. Meanwhile, the leading eigenvalue with the lowest growth rate across the baseline and all control cases also occurs at the optimal angle 
$\psi = 30^\circ$, indicating that reducing the size of the shedding vortices and achieving a balance between the primary and secondary shocks can most effectively suppress the baseline instabilities and stabilize the controlled flow. The resolvent modes obtained from all control cases show that the USL response becomes more unstable and is not sensitive to changes in the angle when 
$\psi <0^\circ$. When 
$\psi > 0^\circ$, the USL becomes less responsive to perturbations, implying the flow is more stable, however, the USL is significantly more sensitive to changes in 
$\psi$. Introducing control with a larger streamwise velocity component is found to amplify the shock and shear layer instabilities. Actuator placement on the bottom surface of the splitter plate constricts the bypass stream and is not as effective as it amplifies the resonant tone and strengthens shock structures. Similar amplifications in the shock strength are observed for the trailing edge surface actuation in cases where 
$\psi < 0^\circ$. This suggests that to achieve success in control efforts, microjet actuation should be introduced upwards into the supersonic main core flow.
The reduction in flow separation, shock strength and the dominant tone throughout the control cases presented in this study reveals the effectiveness of introducing steady blowing at the splitter plate proximity region. This prompts future investigations in achieving similar results through the use of unsteady actuators. Applying the discounted resolvent analysis on the controlled flows can provide insights into how the control configuration can be further optimized through a frequency parameter to achieve the same control outcomes with minimal energy expenditure.
Acknowledgements
We would like to acknowledge the computational resources provided by Syracuse University Research Computing as well as the Syracuse University Fellowship. We also gratefully acknowledge the valuable discussions had with Dr M.N. Glauser.
Funding
This work is supported by the Air Force Office of Scientific Research under award number FA9550-23-1-0019 (Program Officer: Dr G. Abate).
Declaration of interests
The authors report no conflict of interest.
Open access
