2.1. Experimental set-up

The experimental database used in this work comprises of velocity measurements previously reported in Mitchell, Honnery & Soria (Reference Mitchell, Honnery and Soria2013), and acoustic measurements reported in Edgington-Mitchell et al. (Reference Edgington-Mitchell, Honnery and Soria2015).

The elliptical jets studied here are produced from a purely converging nozzle with aspect ratio $AR = 2$. The nozzle is connected to an upstream circular flange, and smoothly transitions from a 70 mm diameter circular section to an ellipse with major and minor axes of 20 and 10 mm, respectively, over a distance of 70 mm, before terminating with a 2.5 mm parallel section. The varying internal length in the azimuthal direction means that the boundary layer will likewise be azimuthally inhomogeneous. The state of the boundary layer at exit has not been established quantitatively, but visualisation of a $M=0.9$ jet from the same nozzle indicates an initially laminar shear layer that rapidly undergoes transition to turbulence. Two jets are considered in this study, at nozzle pressure ratios of $NPR = 2.6$ and $3.4$, which act here as surrogate for the Mach number of these jets – the corresponding ideally expanded Mach numbers for these values of $NPR$ are $M_j \approx 1.25$ and $1.45$.

The velocity measurements are taken over a field of view of $2D_e \times 10D_e$ aligned with the major and minor axes of the ellipse. For each axis and each pressure ratio considered, 1000 statistically independent velocity snapshots were obtained via two-component 2-D particle image velocimetry (PIV). The velocity vectors have a spacing of $0.013D_e$, with a spatial resolution based on minimum interrogation window size of $0.026D_e$, though the through-plane resolution is much lower at $0.1D_e$. A detailed discussion of the PIV measurements is available in Mitchell et al. (Reference Mitchell, Honnery and Soria2013). The velocity fields in the major-axis plane are biased by a slight misalignment of the laser sheet; for the purposes of the analysis in this paper, symmetry is enforced across the centreline to remove this bias.

Acoustic measurements are used in the context of this work only to identify the non-dimensional frequencies associated with the screech tones. A G.R.A.S type 46BE 1/4 in. preamplified freefield microphone located in the minor axis at a polar angle of $\phi = 20^\circ$ from the downstream axis.

The velocity data are further decomposed via a snapshot-based POD (Sirovich Reference Sirovich1987). The autocovariance matrix is constructed from the velocity snapshots $\boldsymbol {V}$ such that

(2.1)\begin{equation} \boldsymbol{R}=\boldsymbol{V}^T\boldsymbol{V}. \end{equation}

The solution of the eigenvalue problem

(2.2)\begin{equation} R\boldsymbol{v}=\lambda \boldsymbol{v}, \end{equation}

yields the eigenpairs ($\lambda,\boldsymbol {v}$) from which the spatial POD modes are constructed as

(2.3)\begin{equation} \boldsymbol{\phi}_n(x,y)=\frac{\boldsymbol{V} \boldsymbol{v}_n(t)}{\|\boldsymbol{Vv}_n(t)\|}. \end{equation}

The coefficients at each time $t$ for each mode $n$ can be expressed as

(2.4)\begin{equation} a_n(t)=\boldsymbol{v}_n(t)\|\boldsymbol{Vv}_n(t)\|. \end{equation}

For a POD mode pair that represents a travelling wave, the two POD modes should be in phase quadrature, which in the absence of noise will produce a perfect circle when their temporal coefficients are plotted against each other.

2.2. Reynolds-averaged Navier–Stokes

The RANS modelling is carried out using the open-source code OpenFOAM within the framework of the density-based High-Speed Aerodynamic (HiSA) (Heyns Reference Heyns2023) solver. Details about the implementation and some verification cases for similar supersonic flow scenarios are given in Macia et al. (Reference Macia, Castilla, Gamez-Montero, Camacho and Codina2019) and Heyns & Oxtoby (Reference Heyns and Oxtoby2014). The AUSM$^+$Up scheme proposed by Liou (Reference Liou2006) is employed for capturing discontinuities in the fluid and a dual time stepping (pseudotime) approach used for iterating the solution towards steady state. Turbulence modelling is performed with the realisable $k$–$\epsilon$ model proposed by Shih et al. (Reference Shih, Liou, Shabbir, Yang and Zhu1995).

The computational domain and mesh with corresponding boundary conditions is shown in figure 1. A fully structured o-grid mesh was devised with $y^+ = 0.1$ at the walls and in order to ensure sufficient grid resolution a grid-convergence study was carried out. At the inlet, the total pressure is ramped up to the desired $NPR$ in order to ensure numerical stability. It is assumed that the operating gas enters the nozzle at ambient temperature (300 K) and the velocity at the inlet is set to Neumann condition (${\rm d}u_i/{\rm d}\kern0.7pt x_i = 0$). At the outlet, standard ambient conditions are set with a fixed pressure of $1 \times 10^5$ Pa and temperature of 300 K.

Figure 1. Schematic of boundary conditions with structured grid for the present simulations.

The ideal-gas model is used for calculating the density of the fluid and the viscosity computed according to Sutherland's law ($\mu _0 = 1.716\times 10^{-5}\ {\rm kg}\ {\rm m}^{-1}\ {\rm s}^{-1}$, $T_0 = 273$ K, $S_\mu = 111$ K).

2.3. Spatial LSA

The linear-stability model used herein is based on work by Lajús et al. (Reference Lajús, Sinha, Cavalieri, Deschamps and Colonius2019) in the study of corrugated jets, and Nogueira & Edgington-Mitchell (Reference Nogueira and Edgington-Mitchell2021) in the study of twin-jets. The coordinate system is given in figure 2. The formulation assumes compressible, inviscid, locally parallel flow, for which the non-dimensional Navier–Stokes equations in cylindrical coordinates are given by

(2.5)$$\begin{gather} \frac{{\rm D}}{{\rm D}t} \left( \frac{1}{\tilde{\rho}} \right) - \left( \frac{1}{\tilde{\rho}} \right)(\boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{u}}) = 0, \end{gather}$$

(2.6)$$\begin{gather}\frac{{\rm D}\tilde{u}_x}{{\rm D}t} + \frac{1}{\tilde{\rho}}\frac{\partial \tilde{p}}{\partial x} = 0, \end{gather}$$

(2.7)$$\begin{gather}\frac{{\rm D}\tilde{u}_r}{{\rm D}t} - \frac{\tilde{u}_\theta^2}{r} + \frac{1}{\tilde{\rho}}\frac{\partial \tilde{p}}{\partial r} = 0, \end{gather}$$

(2.8)$$\begin{gather}\frac{{\rm D}\tilde{u}_\theta}{{\rm D}t} +\frac{\tilde{u}_\theta \tilde{u}_r}{r} + \frac{1}{r} \frac{1}{\tilde{\rho}} \frac{\partial \tilde{p}}{\partial \theta} = 0, \end{gather}$$

(2.9)$$\begin{gather}\frac{{\rm D}\tilde{p}}{{\rm D}t} + \gamma \tilde{p}(\boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{u}}) = 0, \end{gather}$$

with

(2.10)$$\begin{gather} \frac{{\rm D} }{{\rm D}t} = \frac{\partial}{\partial t} + \tilde{u}_x \frac{\partial}{\partial x} + \tilde{u}_r \frac{\partial}{\partial r} + \frac{\tilde{u}_\theta}{r} \frac{\partial}{\partial \theta}, \end{gather}$$

(2.11)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{u}} = \frac{\partial \tilde{u}_x}{\partial x} + \frac{\partial \tilde{u}_r}{\partial r} + \frac{\tilde{u}_r}{r} + \frac{1}{r}\frac{\partial \tilde{u}_\theta}{\partial \theta}, \end{gather}$$

where $\tilde {\boldsymbol {u}}=(\tilde {u}_x,\tilde {u}_r,\tilde {u}_\theta )$ is the velocity vector following figure 2, $\tilde {\rho }$ is the density, $\tilde {p}$ is the pressure and $\gamma$ is the ratio between specific heats. The variables are made non-dimensional by the equivalent diameter $D_{eq}$ (the diameter of a circular nozzle with the same area as the elliptic nozzle as per Hussain & Husain (Reference Hussain and Husain1989)), the free stream density $\rho _\infty$ and the ideally expanded jet velocity $U_j$, which may be obtained from the nozzle-pressure-ratio values using isentropic relations. The equations of motion are then linearised by decomposing variables into mean and perturbation components and neglecting the second-order terms. Azimuthal and radial mean velocities are assumed to be zero in this locally parallel flow assumption. To consider the behaviour of instability waves, a normal-mode ansatz is taken, which assumes that the perturbations take the form of a wave, which assumes the shape $q(r,\theta ) \exp (-\mathrm {i} \omega t +\mathrm {i} \alpha x)$, with $q$ being any of the aforementioned flow variables. In the current spatial stability analysis, $\omega$ is a real number, and $\alpha$ is complex; its real part ($\alpha _r$) corresponds to the wavenumber of the disturbances while its imaginary part ($-\alpha _i$) represents their growth rates. Following Lajús et al. (Reference Lajús, Sinha, Cavalieri, Deschamps and Colonius2019), the resulting set of equations may be rearranged so as to obtain a single equation for the pressure perturbations, given as

(2.12)\begin{equation} \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial p}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2 p}{\partial \theta^2} + A \frac{\partial p}{\partial r} + B \frac{\partial p}{\partial \theta} + C p=0,\end{equation}

with

(2.13)$$\begin{gather} A = \frac{1}{\gamma \bar{p}} \frac{\partial \bar{p}}{\partial r} - \frac{1}{\bar{\rho}} \frac{\partial \bar{\rho}}{\partial r} - \left( \frac{2\alpha}{\alpha \bar{u}_x -\omega} \right) \frac{\partial \bar{u}_x}{\partial r}, \end{gather}$$

(2.14)$$\begin{gather}B = \frac{1}{\gamma \bar{p}} \frac{\partial \bar{p}}{\partial \theta} - \frac{1}{\bar{\rho}} \frac{\partial \bar{\rho}}{\partial \theta} - \left( \frac{2\alpha}{\alpha \bar{u}_x -\omega} \right) \frac{\partial \bar{u}_x}{\partial \theta}, \end{gather}$$

(2.15)$$\begin{gather}C = \frac{\bar{\rho}( \alpha \bar{u}_x -\omega )^2}{\gamma \bar{p}} - \alpha^2. \end{gather}$$

Floquet theory can be applied considering a $2{\rm \pi} /N$ azimuthal periodicity of the mean jet flow (the elliptic jet is ${\rm \pi}$ periodic, thus $N=2$), which allows the pressure-amplitude function to be expressed as

(2.16)\begin{equation} p(r,\theta)=\hat{p}(r,\theta)\,\mathrm{e}^{\mathrm{i}\mu\theta},\end{equation}

where $\mu$ is the Floquet exponent. This allows the domain size to be halved, reducing the computational resources required. The Floquet exponent maps to azimuthal wavenumbers based on the number of lobes $N$ as

(2.17)\begin{equation} m=\{ \ldots, \mu-3N, \mu-2N, \mu-N, \mu, \mu+N, \mu+2N, \mu+3N,\ldots\},\end{equation}

which are all the azimuthal modes $m$ included in the simulations for each value of $\mu$. Substituting equation (2.16) into (2.12) yields

(2.18)\begin{equation} \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial \hat{p}}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2 \hat{p}}{\partial \theta^2} + \hat{A} \frac{\partial \hat{p}}{\partial r} + \hat{B} \frac{\partial \hat{p}}{\partial \theta} + \hat{C} \hat{p}=0,\end{equation}

with

(2.19)$$\begin{gather} \hat{A} = \frac{1}{\gamma \bar{p}} \frac{\partial \bar{p}}{\partial r} - \frac{1}{\bar{\rho}} \frac{\partial \bar{\rho}}{\partial r} - \left( \frac{2\alpha}{\alpha \bar{u}_x -\omega} \right) \frac{\partial \bar{u}_x}{\partial r}, \end{gather}$$

(2.20)$$\begin{gather}\hat{B} = \frac{1}{\gamma \bar{p}} \frac{\partial \bar{p}}{\partial \theta} - \frac{1}{\bar{\rho}} \frac{\partial \bar{\rho}}{\partial \theta} - \left( \frac{2\alpha}{\alpha \bar{u}_x -\omega} \right) \frac{\partial \bar{u}_x}{\partial \theta} + 2\mathrm{i}\mu, \end{gather}$$

(2.21)$$\begin{gather}\hat{C} = \frac{\bar{\rho}( \alpha \bar{u}_x -\omega )^2}{\gamma \bar{p}} - \alpha^2 - \frac{\mu^2}{r^2} + \frac{\mathrm{i} \mu}{r^2} \left[ \frac{1}{\gamma \bar{p}} \frac{\partial \bar{p}}{\partial \theta} - \frac{1}{\bar{\rho}} \frac{\partial \bar{\rho}}{\partial \theta} - \left( \frac{2\alpha}{\alpha \bar{u}_x -\omega} \right) \frac{\partial \bar{u}_x}{\partial \theta} \right]. \end{gather}$$

Figure 2. Coordinate system for LSA formulation.

The equation above may be re-expressed as a polynomial in terms of the wavenumber $\alpha$ as

(2.22)\begin{equation} (\boldsymbol{F}_0 + \alpha \boldsymbol{F}_1 + \alpha^2 \boldsymbol{F}_2 + \alpha^3 \boldsymbol{F}_3) \hat{p} = 0,\end{equation}

where operators $\boldsymbol {F}_n$ are given in Stavropoulos et al. (Reference Stavropoulos, Mancinelli, Jordan, Jaunet, Weightman, Edgington-Mitchell and Nogueira2023). Equation (2.22) can be posed as a nonlinear eigenvalue problem by discretising the domain in both radial and azimuthal directions. This allows the problem to be solved as the matrix equation

(2.23)\begin{equation} \begin{bmatrix} O & I & O \\ O & O & I \\ -\boldsymbol{F}_0 & -\boldsymbol{F}_1 & -\boldsymbol{F}_2 \end{bmatrix} \begin{bmatrix} \hat{p} \\ \alpha \hat{p} \\ \alpha^2 \hat{p} \end{bmatrix} = \alpha \begin{bmatrix} I & O & O \\ O & I & O \\ O & O & \boldsymbol{F}_3 \end{bmatrix} \begin{bmatrix} \hat{p} \\ \alpha \hat{p} \\ \alpha^2 \hat{p} \end{bmatrix}, \end{equation}

where $I$ and $O$ are identity and null matrices and $\boldsymbol {F}_0$, $\boldsymbol {F}_1$, $\boldsymbol {F}_2$ and $\boldsymbol {F}_3$ are the matrices representing the discretised operators in (2.22). This can be solved numerically for the eigenvalue $\alpha$ and eigenfunction $\hat {p}$, which can then be used to calculate the pressure-amplitude function using (2.16). Once the pressure eigenfunction is determined, the linearised governing equations can be used to determine the velocity and density eigenfunctions (Lajús et al. Reference Lajús, Sinha, Cavalieri, Deschamps and Colonius2019).

A Dirichlet boundary condition setting the pressure-amplitude function to zero in the far-field is implemented, since pressure fluctuations associated with instability waves should be vanishingly small far away from the jet. However, since the domain is finite, the pressure eigenfunction is set to zero at the maximum radial coordinate $r_{max}$, as to approximate the exponentially decaying behaviour of that quantity. At the centre of the jet, boundary conditions are set as (Gudmundsson Reference Gudmundsson2009; Nogueira & Edgington-Mitchell Reference Nogueira and Edgington-Mitchell2021)

(2.24)\begin{equation} \left[ r^2\frac{\partial^2 p}{\partial r^2} + r\frac{\partial p}{\partial r} - \mu^2 p \right]_{r\to0} = 0. \end{equation}

A cylindrical-polar mesh was used in these calculations, which consists of a radial-mesh component and an azimuthal-mesh component. Azimuthal points are equispaced, and azimuthal derivatives were computed using Fourier differences (Weideman & Reddy Reference Weideman and Reddy2000), whereas the radial mesh is constructed from Chebyshev polynomials, which are then used to compute radial derivatives (Trefethen Reference Trefethen2000). A transformation from Bayliss & Turkel (Reference Bayliss and Turkel1992) was applied to concentrate the radial points around the shear layer. Both the validation of the method and a convergence study are available in the Appendix (A).

Linear stability analysis aims to evaluate particular solutions for the linearised Navier–Stokes equations under a locally parallel assumption by neglecting the nonlinear effects of turbulence on the wave-like disturbances. The method is also a function of the choice of base-flow around which one linearises the equations; in the present case, these base flows are given either by the experimental mean flow obtained by limited PIV data, or by an approximate solution for the time-averaged flow under a turbulence closure model given by RANS. The application of linear-stability approaches is formally justified in the context of resolvent analysis, as described in Beneddine et al. (Reference Beneddine, Sipp, Arnault, Dandois and Lesshafft2016); the authors showed that the analysis generally produces faithful results for convectively unstable, slowly varying flows such as cold jets. Even though shock-containing jets are not slowly varying, application of linear stability in this case still leads to a good description of the waves involved in resonance (see, for instance, Edgington-Mitchell et al. (Reference Edgington-Mitchell, Jaunet, Jordan, Towne, Soria and Honnery2018), Gojon et al. (Reference Gojon, Bogey and Mihaescu2018); Edgington-Mitchell et al. (Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021) and Nogueira & Edgington-Mitchell (Reference Nogueira and Edgington-Mitchell2021)).

3.1. Experimentally derived two-dimensional base flow and fitted profiles

Unlike axisymmetric jets, the overall mean flow characteristics of elliptical jets cannot be educed by velocity measurements in a single PIV plane; the tri-dimensional nature of the elliptical jet dictates a fully three-dimensional (3-D) measurement of the velocity field (such as that achieved for a rectangular jet by Sellappan & Alvi (Reference Sellappan and Alvi2022)) is necessary to characterise the shear layer development in both $x$ and azimuth. In the absence of such a measurement, the velocities in the minor and major axes of the jet were measured, which are shown in figure 3. In both $NPR$ cases, complex shock-cell structures are observed, especially in the major axis; a low-velocity region is present in that axis (at $z/D_{eq}=0$), which seems to be necessary to match the centreline velocity in both axes. The jet expands more strongly in the minor axis, but the same periodicity is observed in both axes. While shear layer and potential core of the higher $NPR$ jet develops more slowly, the $NPR=2.6$ jet experiences a somewhat stronger variation around $x/D=5$, where the core of the jet shrinks substantially in the major axis, and the shear-layer expands more quickly in the minor. These features seem to indicate axis switching at that position for the lower-$NPR$ jet. It is possible that axis switching also occur for the higher-velocity case, but farther downstream from the measurement window.

Figure 3. Mean streamwise velocity fields obtained from PIV for $NPR = 2.6$ (a,c) and $NPR = 3.4$ (b,d). Data taken from Mitchell et al. (Reference Mitchell, Honnery and Soria2013).

For a given base flow, LSA can calculate the two-dimensional (2-D) velocity, pressure and density eigenfunctions, with their corresponding growth rate and wavenumber. Therefore, if the base flow of the experimental PIV data can be determined, the velocity eigenfunctions on the minor axis from LSA can be compared with the experimentally determined coherent velocity fluctuations ($\psi$) from POD (Nogueira & Edgington-Mitchell Reference Nogueira and Edgington-Mitchell2021). In addition to the base flow, LSA requires inputs of frequency (in its non-dimensional form, Strouhal number), Mach number $M$ and static-temperature ratio (or the jet-to-ambient temperature ratio). Isentropic relations were used to calculate Mach number and static-temperature ratio for each $NPR$ (the experimental operating condition). To apply LSA to the experimentally derived data requires an additional modelling step, however. The velocity profile is only known along the major and minor axes, not for the full 2-D $(r,\theta )$ velocity field. Therefore, for a given axial location, the base flow must be fitted to some sort of analytical mean-flow model that captures the characteristics of the jet in both major and minor axes, modelling the region in between. An expression based on Morris & Bhat (Reference Morris and Bhat1995) was used, which assumes an elliptical coordinate system as

(3.1)$$\begin{gather} y= a_e \cosh\varrho \cos \phi, \end{gather}$$

(3.2)$$\begin{gather}z= a_e \sinh\varrho \sin \phi, \end{gather}$$

where $a_e$ is the location of the foci of the ellipse (at $\pm a_e$), and $\varrho$, $\phi$ are the elliptic radial and azimuthal coordinates, respectively. An elliptical base flow similar to the one proposed by Morris & Bhat (Reference Morris and Bhat1995) is defined as

(3.3)\begin{equation} \frac{U}{U_j} = \begin{cases} U_c, & \text{if } \eta < 0, \\ U_c \, \mathrm{e}^{- \eta^2}, & \text{if } \eta \geq 0,\\ \end{cases}\end{equation}

with

(3.4)$$\begin{gather} \eta=f(\phi) \left( \frac{\text{asin}\varrho - h_B}{b_B} \right), \end{gather}$$

(3.5)$$\begin{gather}f(\phi) = \left( \frac{1+f_0}{2} \right) - \left( \frac{1-f_0}{2} \right) \cos(2\phi), \end{gather}$$

(3.6)$$\begin{gather}f_0=\frac{b_B}{\sqrt{(h_A+b_A)^2-a_e^2}-h_B}. \end{gather}$$

In the equations above, $h_A$ and $h_B$ are the jet core width of the major and minor axes, respectively, and $b_A$ and $b_B$ are the half-widths of the shear layer on the major and minor axes, respectively (the distance at which the velocity decreases to half the centreline velocity). The parameter $U_c$ represents the centreline velocity, which was obtained by smoothing out the shock-cells from the data at $y=z=0$. The aspect ratio of the flow is defined as

(3.7)\begin{equation} AR=\frac{h_A+b_A}{h_B+b_B}. \end{equation}

For each axial station of the PIV data, (3.3) was fitted to the data using a nonlinear least squares method, from which the potential core width $(h_A, h_B)$ and the shear-layer half-widths $(b_A, b_B)$ on the major and minor axes, respectively, were determined, allowing a 2-D velocity field to be generated. The fitted potential-core widths, half-widths and the aspect ratio of the PIV flow fields are given in figure 4. As expected from figure 3, the potential-core width decreases with streamwise distance for both cases, while the shear-layer half-width increases. Close to the axis switching point, the core width in the major axis decreases more strongly, highlighting the changes in the mean flow around that position. More interestingly, the aspect ratio of the jets does not start exactly at $2$; at $x/D_{eq}=0.5$, its value is around $1.7$ for both jets. This discrepancy with the nominal $AR$ of the jet may be due to the definition of a flow aspect ratio, which is somewhat arbitrary. However, a clear trend of reduction is observed as $x/D_{eq}$ is increased, highlighting how the jet becomes increasingly circular as the shear layer develops (at least if one considers the mean flow to be well described by the equation of the ellipse as the jet evolves). The axis switching point ($AR=1$) is around the position where the mean flow is observed to change the most for $NPR=2.6$, as mentioned previously; the axis switching point is located farther downstream for $NPR=3.4$, being not as noticeable as in the lower-speed case.

Figure 4. Fitted core widths (a,d) and shear layer half-widths (b,e) to PIV profiles, and the resulting aspect ratio (c,f) of the flow for $NPR=2.6$ (a–c) and $3.4$ (d–f).

A comparison of the PIV velocity profile and the ideally expanded fitted profile on the major and minor axes are shown in figure 5, where the relevant first derivatives of the mean flow are also shown. It is clear that, while the original velocity profiles display complex features in the potential core resulting from the presence of shock-cells, the fitted profile seems to capture the overall shape of the flow, especially outside of the core. This is confirmed by the gradients shown in figure 5, which show good agreement between the fit and the experiments. A similar level of agreement is observed in both major and minor axes, at upstream and downstream regions. As shocks are generally undesirable in the LSA computations, this profile choice seems to be well suited for the analysis at hand.

Figure 5. Comparison between fitted velocity profiles and actual velocity profiles and their gradients at $x/D_{eq}=0.4$ and $x/D_{eq}=4.0$ for $NPR = 2.6$ (a,c) and $3.4$ (b,d).

3.2. The RANS-derived 3-D base flow

The objective of using RANS in this work is two-fold: first, it provides the actual streamwise velocity distribution in $y$–$z$, which is lacking in the current PIV set-up; this allows us to check if the current fit choice is appropriate to represent the development of the jet in the downstream direction. A secondary objective is to provide data for a sensitivity analysis. As the RANS results capture some of the features of the experimental data, one may use these results to evaluate what are the essential features of the mean flow that lead to good predictions of the waves responsible for screech.

A comparison between the streamwise velocity at the centreline for both $NPR$ is shown in figure 6. Clearly, the simulation results do not capture all the features from the experimental data, with the most striking difference being the rapid decay of the centreline velocity in the PIV for $NPR=2.6$, which is absent in the RANS. While the disparity farther downstream is quite severe, the key LSA calculations for this pressure ratio presented in later sections are confined to the region $0 \leq x/D_e \leq 4$, where the data agree quite well. The rapid decay of shock cells at $NPR = 2.6$ is likely driven by a violent flapping motion associated with a strong screech tone at this condition (Edgington-Mitchell et al. Reference Edgington-Mitchell, Honnery and Soria2015). The change in spacing of the shock for $NPR=3.4$ is also misrepresented in the simulation, even though the decay of the shocks follows the same overall behaviour. Despite these discrepancies, the first shock-cell wavelength and the maximum amplitude of the shocks are well represented by the simulation for both $NPR$. For the $NPR=2.6$ case, the overall decay of the shocks and the decrease in spacing farther downstream is also observed in both experiments and simulation, even though the minima of velocity induced by the shocks are not well captured. Overall, it is clear that the simulation has only captured the flow in a qualitative sense; however, it is well known that results from RANS are dependent on the choice of turbulence model and their underlying parameters, and no attempt has been made to tune these parameters here. As the simulation cannot capture the effects of screech, which include a significant amplification of shear-layer growth, it is unsurprisingly than an untuned turbulence model produces only qualitative agreement. However, the purpose of the RANS here is not to exactly reproduce a ‘real’ flow, but to produce internally consistent solutions that can act as an input to the LSA calculations, as an alternative to the empirical mean-flow fit used with PIV data. Thus, the current simulation results may only be seen as an approximation of the real data, but one that may give insight about the blind spots induced by the current experimental set-up, namely the absence of cross-plane data.

Figure 6. Comparison between streamwise velocity magnitudes in the centreline from PIV and RANS for $NPR=2.6$ (a) and $NPR=3.4$ (b).

The fitted mean flow profile and the RANS are plotted for $x/D_{eq}=0.5$, $2$ and $4$ for both $NPR$ in figures 7 and 8. Very close to the nozzle, the mean flows are remarkably similar for both cases, with the exception of the shock-cell signature in the core of the jet visible in the RANS result. The level of agreement deteriorates with increasing streamwise position for the $NPR=3.4$ case, especially at $x/D_{eq}=4$, where the mean flow does not resemble a perfect ellipse anymore; at this position, the RANS suggests that the real mean flow will have a more intricate shape. For the lower $NPR$ case, there is less deviation from an elliptical shape, but there is significantly greater disparity in shear layer thickness. Consistent with the previously described mismatch in centreline decay, the shear layer of the fitted profile is also significantly thicker than the simulation. This is once again likely due to the violent flapping motion induced by screech not being captured by the RANS. Despite these differences in shear layer thickness, the qualitative match of shape is excellent for the streamwise positions considered in figure 7, supporting the choice of velocity profile used in the fit.

Figure 7. Streamwise velocity profiles at $x/D_{eq}=0.5$ (a,d), $2$ (b,e) and $4$ (c,f) for $NPR = 2.6$. Both experimentally fitted (a–c) and RANS (d–f) profiles are shown.

Figure 8. Streamwise velocity profiles at $x/D_{eq}=0.5$ (a,d), $2$ (b,e) and $4$ (c,f) for $NPR = 3.4$. Both experimentally fitted (a–c) and RANS (d–f) profiles are shown.

3.3. Wave-like structures in the experimental data

Fluctuating-velocity fields were determined by subtracting the mean flow from each snapshot obtained from PIV. Proper orthogonal decomposition was then performed on the fluctuating-transverse-velocity fields on the major and minor axes (as described in § 2.1), yielding the spatial shape of the most-energetic coherent structures in the flow and their relative energies. The transverse component was chosen as to avoid the influence of shear-thickness modes (Weightman et al. Reference Weightman, Amili, Honnery, Soria and Edgington-Mitchell2018; Nogueira & Edgington-Mitchell Reference Nogueira and Edgington-Mitchell2021). The relative energy of the POD modes on the major and minor axes for each $NPR$ is presented in figure 9.

Figure 9. Transverse velocity fluctuation POD mode energies for $NPR = 2.6$ (a) and $NPR = 3.4$ (b).

It is evident that for both major and minor axes the first two modes contain significantly more energy than all other modes, with the exception of the major axis for $NPR = 2.6$ where the first four modes are of similar energy. On the minor axis, the first two modes account for $35\,\%$ of total energy for $NPR = 3.4$, and $23.6\,\%$ for $NPR = 2.6$, whereas on the major axis the first two modes comprise $10.3\,\%$ of total energy for $NPR=3.4$ and $6.8\,\%$ for $NPR=2.6$. Oberleithner et al. (Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011) suggests that if the energy of the first two POD modes are substantially higher than the other modes, and they are of similar magnitude, the modes may represent a travelling wave. As described in § 2.1, for a pair of modes in phase quadrature, the phase portrait of their temporal coefficients should form a circle. This is shown in figure 10, which demonstrates that the phase portrait of the first and second minor axis modes for each $NPR$ forms a circle, whereas the phase portrait of the first and second modes on the major axis does not. This further suggests that the first two POD modes of the minor axis may comprise of a single wave, whereas wave-like behaviour is less pronounced on the major axis. This is consistent with the tendency of jets to ‘flap’ about their minor axis.

Figure 10. Phase portraits of first two POD modes at each axis for $NPR = 2.6$ (a) and $NPR = 3.4$ (b).

Figures 10 and 11 together confirm that both jets are resonating in a flapping mode in the minor axis, leading to symmetric velocity disturbances in the transverse direction. In order to improve mode convergence and facilitate comparison with modes predicted by LSA, the POD was recomputed using symmetry-imposed fields, as in Nogueira & Edgington-Mitchell (Reference Nogueira and Edgington-Mitchell2021). Following the results of figure 11, only even modes in transverse velocity will be analysed. The imposition of symmetry does not significantly change the mode shapes or energies, but does slightly improve convergence.

Figure 11. Transverse velocity fluctuations associated with minor axis POD modes 1 (a,b) and 2 (c,d) for $NPR=2.6$ (a,c) and $NPR=3.4$ (b,d).

Since snapshot POD gives no time resolution, a complex travelling wave can be constructed as (Oberleithner et al. Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011; Edgington-Mitchell et al. Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021)

(3.8)\begin{equation} u_c(x,y,t)=\hat{a}[ \phi_1(x,y)+\mathrm{i} \phi_2(x,y) ]\, \mathrm{e}^{-\mathrm{i}\omega t}, \end{equation}

where $\hat {a} \,\mathrm {e}^{-\mathrm {i}\omega t} =a_1-\mathrm {i} a_2$. The equation above provides the spatial support of the most amplified travelling wave in the flow, which is expected to be associated with the screech phenomenon – this wave-like mode is now represented as a complex mode. The absolute value of the resulting structure is shown in figure 12. These plots show that the most-energetic coherent structure is dominated by a standing wave, which is especially prevalent for $NPR = 3.4$. The standing wave is a signature of the jet undergoing screech, being caused by a superposition between the upstream travelling GJM and downstream travelling KH instability (Panda Reference Panda1999; Edgington-Mitchell et al. Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021). The standing-wave is observed across a large region of the space, which includes the centreline of the jet, as in the flapping/helical mode studied in Nogueira et al. (Reference Nogueira, Jordan, Jaunet, Cavalieri, Towne and Edgington-Mitchell2022b).

Figure 12. Absolute value of the resonance mode reconstructed using POD modes 1 and 2 on minor axis for $NPR = 2.6$ (a) and $NPR = 3.4$ (b). Only the transverse velocity is used in the reconstruction.

A spatial Fourier transform can be performed to deduce the most energetic wavenumbers $k_x$ of the waves that produce the standing wave. The wavenumber spectrum of the Fourier-transformed coherent-velocity-amplitude function is presented in figure 13. As in previous works (Edgington-Mitchell et al. Reference Edgington-Mitchell, Jaunet, Jordan, Towne, Soria and Honnery2018, Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021; Nogueira & Edgington-Mitchell Reference Nogueira and Edgington-Mitchell2021), phase and group velocity are assumed to have the same sign. Therefore, for a positive frequency, positive wavenumbers correspond to downstream-travelling waves and negative wavenumbers correspond to upstream-travelling waves. Figure 13 demonstrates two prominent peaks on each wavenumber spectra; a peak for $k_x>0$ (downstream travelling) and a peak for $k_x<0$ (upstream travelling). The upstream-travelling peak is at a wavenumber corresponding to a phase velocity slightly slower than the speed of sound. This is consistent with the current understanding of screeching jets, where the downstream travelling wave is the KH instability, whereas the upstream travelling wave is the GJM (Edgington-Mitchell et al. Reference Edgington-Mitchell, Li, Liu, He, Wong, Mackenzie and Nogueira2022). Table 1 shows the peak wavenumber for each of the waves in figure 13 and the range of wavenumbers that comprise each peak, computed as the wavenumbers where the energy decays to $10\,\%$ of the peak. A secondary downstream peak is also visible in figure 13 very close to the dominant KH peak. The origin of this peak is unclear, but it is possibly related to the limitations of POD in isolating the structures at a single frequency, or to a secondary structures associated with the development of the vortices. The lack of time-resolved data hinders a clearer identification of the physical mechanism behind this peak in the present case, and further analysis of the possible causes will be left as a future work, as this secondary peak does not affect our current interpretation of the phenomenon.

Figure 13. Wavenumber spectra of the dominant travelling wave for $NPR = 2.6$ (a) and $NPR = 3.4$ (b).

Table 1. Wavenumber range of resulting KH and GJM extracted from POD.

To isolate the structure of one of each of these waves, all velocity amplitude functions outside of the wavenumber range of each wave (provided in table 1) are set to zero. Then an inverse spatial Fourier transform can be performed to determine the spatial structure of the waves represented in the wavenumber spectrum. The resulting structure of the GJM and the KH instability wave are presented in figure 14. It is evident that, when comparing each $NPR$, the equivalent KH and the GJMs have the same overall structure, although the peak magnitude for each of these waves occurs farther downstream for the $NPR=3.4$ jet. This is consistent with figure 12, which shows that the standing wave pattern for the $NPR=3.4$ jet reaches a higher amplitude farther downstream than the $NPR=2.6$ jet. Interestingly, both waves follow the same behaviour as the flapping/helical mode in circular jets shown in Edgington-Mitchell et al. (Reference Edgington-Mitchell, Wang, Nogueira, Schmidt, Jaunet, Duke, Jordan and Towne2021), at least in this plane of measurement. Particularly, both waves have peaks at the centreline of the jets. For the KH wavepacket, it is expected that the disturbances peak at the shear layer if it is sufficiently thin; however, as this mode develops downstream, the peak in the transverse velocity component for flapping-like disturbances (with dominant azimuthal wavenumber $|m|=1$) may shift to regions closer to the centreline, as will be shown shortly using LSA. One should note that the mode reconstruction in figure 14 is performed by considering waves in the wavenumber range shown in table 1 (computed as the wavenumbers where the energy decays to $10\,\%$ of the peak). As expected, the spatial shapes of the modes are somewhat dependent on that threshold value, but they remain qualitatively the same. The streamwise positions of the peaks of both waves change less than $5\,\%$ when the threshold is varied between $10\,\%$ and $50\,\%$, confirming the robustness of the approach.

Figure 14. Filtered GJM (a,c) and KH instability (b,d) from reconstructed POD wave for $NPR = 2.6$ (a,b) and $NPR = 3.4$ (c,d). Only the transverse velocity component is shown.

4.1. The KH mode

As stated in § 3.3, coherent velocity fluctuation associated with travelling waves was only observed about the major axis (along the minor axis). Morris & Bhat (Reference Morris and Bhat1995) have shown using LSA that the KH instability splits into two different Mathieu function modes, the ‘cosine-elliptic’ ($ce_n$) and ‘sine-elliptic’ ($se_n$). It is well known that the $se_1$ KH pressure fluctuation mode has odd symmetry about the major axis, and will cause ‘flapping’ velocity fluctuations along the minor axis, as shown in Morris & Bhat (Reference Morris and Bhat1995) and in the Appendix (A). Since screeching elliptical jets are known to flap along that axis, the $se_1$ mode is likely to be the mode that is observed in the experiments; therefore only the $se_1$ KH mode will be analysed using the LSA. One should note that the flapping mode may not be the most unstable mode at a given frequency and streamwise station, but it is the mode involved in resonance in the present case, which justifies focusing on this specific mode in the LSA. The frequency chosen in the analysis is also the screech frequency for each case as determined via acoustic measurements. Comparison of the radial velocity amplitude function from LSA and the POD at various axial locations, at the screech frequency for each $NPR$ is shown in figure 15 for $NPR = 2.6$ and $NPR = 3.4$, using both the PIV-based mean flow and the one obtained from RANS.

Figure 15. Comparison between the absolute value of transverse velocity component obtained from the filtered POD and LSA for the KH mode at several streamwise stations and at the screech frequency. Results were generated at $x/D_{eq}=0.5$, $1.5$ and $2.5$ for $NPR = 2.6$ (a–c) and at $x/D_{eq}=1$, $2$ and $3$ for $NPR=3.4$ (d–f).

Figure 15 demonstrates that the distribution of radial velocity fluctuation along the minor axis is similar between LSA and POD at nearly all axial locations when analysing the PIV data. One should note that the streamwise positions were chosen as to include positions where the KH mode has high amplitudes in the POD, while still being unstable in the Navier–Stokes spectrum. In figure 15, both the radial locations and relative magnitudes of the maxima are close to the experiments whereas the minima of the LSA tend to be slightly larger than that of the POD. Though the agreement is generally good, it is expectedly imperfect; the match at $x/D_{eq}=2.5$ for $NPR=2.6$ for instance is relatively poor. Results using the RANS profile agree with the PIV ones close to the nozzle, but diverge from both the POD and the PIV-based ones farther downstream (while still keeping the same overall features). This may be due to the presence of the shocks, which are known to affect the characteristics of the predicted modes (Nogueira et al. Reference Nogueira, Jaunet, Mancinelli, Jordan and Edgington-Mitchell2022a), or due to the differences in shear-layer thickness. Despite some discrepancies, these results are still encouraging, as the general features of the POD mode can still be predicted by the linear model for two different choices of mean flow. Given one mean flow is based on a fit that assumes an ellipse, while the other includes 3-D effects and shocks, the qualitative similarity in the predicted mode shapes suggests that details of the mean streamwise velocity may have only second-order effects on the KH mode dynamics.

Figures 16 and 17 illustrate the wavenumber and growth rates of the KH instability predicted by the LSA at various downstream stations using both PIV and RANS mean flows. For each $NPR$ in figure 16, the growth rate decreases in the downstream direction. This is consistent with the findings of Nogueira et al. (Reference Nogueira, Jaunet, Mancinelli, Jordan and Edgington-Mitchell2022a) for a circular jet. This is likely due to the thickening of the shear layer (see figure 3), which smooths out velocity gradients, thus stabilising the KH mode, as shown in a range of previous works (see, for instance, Michalke (Reference Michalke1984) and Morris (Reference Morris2010)). In figure 16, the farthest downstream point provided for the KH instability at each $NPR$ is the last axial location for which the mode was found to be unstable. The process of obtaining growth rates from the PIV data is less straightforward and therefore no comparison is drawn. Growth rates of the KH mode have also been shown to be less important in defining the frequency of the resonance process than the actual wavenumber of this mode; they do, however, impact the spatiotemporal growth rates of the resonant mode (Nogueira et al. Reference Nogueira, Jordan, Jaunet, Cavalieri, Towne and Edgington-Mitchell2022b). While not directly related to the resonance model used in the present work, the long-range resonance model of Mancinelli et al. (Reference Mancinelli, Jaunet, Jordan and Towne2021) uses the growth rates of the KH to define positions where the resonance might be sustained in space, which allows for the elimination of spurious-tone predictions in that model.

Figure 16. Growth rate of KH mode for $NPR = 2.6$ (a) and $NPR = 3.4$ (b) as function of streamwise station. Quantities computed at the screech frequency as function of streamwise distance from the nozzle.

Figure 17. Wavenumber of KH mode for $NPR = 2.6$ (a) and $NPR = 3.4$ (b) as function of streamwise station. Quantities computed at the screech frequency as function of streamwise distance from the nozzle.

Similarly, for each $NPR$ in figure 17 it is evident that the wavenumber increases before reaching a maximum, before decreasing. Phase velocity is inversely proportional to wavenumber and therefore this trend implies an initial decrease in phase velocity before reaching a minima, followed by an increase in phase velocity, consistent with the measurements of Mercier, Castelain & Bailly (Reference Mercier, Castelain and Bailly2017). The peak wavenumber reached in the $NPR = 2.6$ jet is observed to occur closer to the nozzle than in the $NPR = 3.4$ jet; this could be associated with the fact that maximum wave amplitude is also reached closer to the nozzle in the $NPR = 2.6$ case. In the POD, the peak wavenumber for $NPR=3.4$ is $2.6$, whereas for $NPR=2.6$ the peak wavenumber is $3.7$. These values are represented by the red lines in figure 17, whilst the horizontal black lines show the wavenumbers for which a $90\,\%$ reduction in the peak-velocity amplitude was observed in the spatial spectrum of the screech mode (see figure 13 and table 1). It is evident that, for each $NPR$, all wavenumbers determined from the LSA fall within the same range determined from the POD, and they are remarkably close to the peak wavenumber for both Mach numbers. The similarity between the structure of the observed transverse velocity fluctuation waves in the POD and the LSA and the quantitative agreement between their wavenumbers, suggest that LSA accurately predicts the behaviour of the $se_1$ KH instability.

The wavenumbers determined by analysis of the RANS data and the PIV data are extremely similar, although the change in wavenumber with axial direction has fluctuations, likely due to the more pronounced effects of the shocks in the RANS solution. The growth rates determined by the analysis of the RANS data are very similar to the PIV data closer to the nozzle (from $x/D = 0$ to $x/D = 1.5$) for $NPR=2.6$, whereas farther downstream the RANS data produces a higher growth rate. The difference is likely due to the fact that the velocity profiles of the PIV and RANS data match quite well when close to the nozzle, and become less similar downstream as the flow profile becomes less elliptic; as shown in figure 7, the shear layer in the PIV thickens more rapidly than in the RANS. These differences are less pronounced in the $NPR=3.4$ case, given the stronger similarities between the mean flows used in this case.

4.2. The GJM

The GJM is known to be the upstream-travelling wave that closes the resonance loop in jet screech. As such, the upstream-travelling wave observed in the POD is expected to take the form of a GJM, and the ability of the LSA to predict the characteristics of this wave will be assessed. Tam & Hu (Reference Tam and Hu1989) (see also Bogey (Reference Bogey2021)) have previously shown, using a vortex sheet model, that the GJM can exist for various azimuthal wavenumbers in an axisymmetric jet. For an elliptical jet, the GJM is expected to follow a similar ‘cosine-elliptic’ and ‘sine-elliptic’ identified in the KH mode, as highlighted in the previous section. The pressure eigenfunctions of GJMs shown in figure 18 below were determined by performing LSA on the PIV data.

Figure 18. Absolute value of pressure eigenfunctions of the GJM for $NPR=2.6$ and $x/D_{eq}=2$. Modes computed for $St=0.2$ (a), $St=0.36$ (b), $St=0.42$ (c).

Figure 18 shows a single $ce_0$ and two flapping modes associated with $ce_1$ and $se_1$ Mathieu functions; note that these modes were computed at different frequencies, as they only exist for a limited range of frequencies and wavenumbers in the spectrum (Tam & Hu Reference Tam and Hu1989). The structures in figure 18 resemble the ones computed for circular (Nogueira et al. Reference Nogueira, Jordan, Jaunet, Cavalieri, Towne and Edgington-Mitchell2022b) and twin (Nogueira & Edgington-Mitchell Reference Nogueira and Edgington-Mitchell2021) jets, with the main difference being the fact that flapping modes have now a clear elliptical shape. In the POD data, wave-like behaviour was only observed in the velocity fluctuations about the major axis, and therefore the GJM observed in the POD is likely of $se_1$ symmetry, as was the KH instability; in fact, it has been observed that this symmetry is usually active around the screech frequency. As such the remainder of the analysis of the GJM will focus on the $se_1$ symmetry.

The transverse velocity on the minor axis was calculated for the $se_1$ GJM from the pressure eigenfunction for several axial stations close to the peak of the POD, and then compared with the experiments. One should note that, given the small range of frequencies for which this mode is neutrally stable, it may not exist for both PIV and RANS mean flows at the same frequency and streamwise position. The comparisons are shown in figure 19, though as mentioned not every panel has predictions from both PIV and RANS. It is evident that, regardless of the $NPR$, axial location or data analysed (PIV or RANS), the structure of the $se_1$ GJM transverse velocity component on the minor axis is somewhat invariant and consists of a large, sharp peak in the centre, and a smaller peak around the location of the shear layer. For each axial location and $NPR$, the centre peak of the mode in both the LSA and the POD match almost perfectly, whereas the relative magnitude and the width of the secondary peak changes slightly with downstream direction.

Figure 19. Comparison between the absolute value of transverse velocity component obtained from the filtered POD and LSA for the GJM at several streamwise stations and at the screech frequency. Results were generated at $x/D_{eq}=3$, $3.5$ and $4$ for $NPR = 2.6$ (a–c) and at $x/D_{eq}=4.5$, $5.5$ and $6.5$ for $NPR=3.4$ (d–f).