Modal analysis of screeching elliptical jets

Abstract

This work presents models for the behaviour of both upstream- and downstream-travelling waves in screeching elliptical jets. Proper orthogonal decomposition is performed on experimental velocity data in both the major and minor axis planes, for an aspect ratio $AR=2$ converging elliptical jet operating at nozzle pressure ratios of $2.6$ and $3.4$. From this decomposition, the radial and axial structure of the guided-jet mode (GJM) and the Kelvin–Helmholtz instability are educed. Linear-stability analysis (LSA) is performed using both the experimentally obtained mean flow, and one obtained using Reynolds-averaged Navier–Stokes (RANS) at matched conditions. It is shown that the wavenumber predicted by LSA for both waves are within the range of experimentally observed wavenumbers. Furthermore, the model accurately predicts the structure of these waves at multiple axial locations, using either the experimental or RANS mean flow. Most critically, it is demonstrated that the GJM is only predicted to be neutrally stable at the screech frequency for a relatively limited streamwise domain, the size and location of which is dependent on the nozzle pressure ratio. A comparison with the amplitude envelope for the GJM extracted from the experimental measurements indicates that the maximum fluctuations associated with the GJM are collocated with this region of the flow that is predicted to support the GJM. While there have been extensive discussions about the frequency dependence of the GJM, this is the first demonstration that its existence is highly dependent on streamwise position within the flow.

1. Introduction

Jet noise is a significant problem, as it can cause permanent hearing loss, structural damage and disruption to those who live near flight paths and airports. As such, increasingly stringent noise regulations are being enforced, thus increasing the demand for quieter jet engines. Lighthill's seminal work (Lighthill 1952) determined that acoustic power output for a jet is proportional to the exhaust velocity to the eighth power. Thus, jet noise reductions have been achieved by implementing ever-increasing high-bypass ratio jet engines, which have higher mass flow rates, and lower exhaust velocities. However, the future viability of increasing the bypass ratio to decrease jet noise is limited. Increasing the bypass ratio requires larger engines to maintain thrust, and engines on commercial jets cannot get much larger without redesigning the aircraft. Creative noise-reduction techniques are now required to achieve further noise improvements.

Orderly wave-like structures have been attributed to noise generation in jets since the 1960s. Some of these waves amplify over some distance or time as they travel downstream, thus being referred to as instability waves. Instability waves in axisymmetric jet flows were first studied using temporal linear-stability analysis (LSA) by Batchelor & Gill (1962) to determine the frequencies and the temporal growth rates of instability waves that are supported by the flow. The relation of coherent structures to jet noise, first postulated by Mollo-Christensen (1967), was further investigated by Crow & Champagne (1971) who compared the instability waves predicted using temporal LSA with those observed in experiments. Spatial LSA (which focuses on analysing structures growing in space, instead of time) has successfully been used to predict the characteristics of instability waves in finite-thickness jets; Michalke (1971) demonstrated close agreement with the experimental phase-velocity measurements of Crow & Champagne (1971). Since then, LSA has been used consistently to predict a range of characteristics of the coherent structures associated with jet noise, being used extensively in noise prediction models (Jordan & Colonius 2013; Cavalieri, Jordan & Lesshafft 2019).

Non-axisymmetric nozzles have the potential to reduce jet noise due to their enhanced mixing characteristics (Chauhan, Aravindh Kumar & Rathakrishnan 2015; Chauhan, Kumar & Rathakrishnan 2016) or due to their ability to disrupt the most unstable modes of instability waves in the flow (Crighton 1973). The closest non-axisymmetric nozzle geometry related to circular nozzles is elliptical nozzles. These nozzles have been observed to radiate less noise than equivalent circular jets at the same operating condition (Kinzie & McLaughlin 1999), confirming Crighton's hypothesis. However, the specific mechanism in which this reduction is achieved has not been exposed, and LSA applied to elliptical jets may allow the underlying physics to be better understood. It is also unclear how universal these reductions are, with studies on subsonic elliptical jets (Amaral et al. 2023), and elliptical nozzles (Viswanathan 2023) not demonstrating clear noise reductions.

Temporal and spatial LSA was first extended to elliptical jet flows in the limit of infinite aspect ratio (Crighton 1973) ascertaining that ‘flapping’ spatial instabilities are dominant. Spatial LSA was extended to finite aspect ratio ($AR$) elliptical jets (Morris 1984, 1988; Morris & Bhat 1995) determining that the growth of the equivalent axisymmetric ‘varicose’ mode in circular jets was damped by the elliptical geometry, although ‘flapping’ modes were observed to become more unstable. Since the paper by Morris & Bhat (1995), few publications have been made extending the use of LSA on elliptical jets (a noteworthy exception being Suzuki, Nogueira & Edgington-Mitchell (2023)). Although Morris was able to explore how Mach number, shear layer thickness and aspect ratio affects the characteristics of instability waves, there was no comprehensive comparison drawn between the LSA model and experimental data. In depth comparisons between experimental data using proper orthogonal decomposition (POD) and LSA have been performed on swirling jets by Oberleithner et al. (2011) and on twin jets by Nogueira & Edgington-Mitchell (2021) which has yet to be extended to elliptical jets. In addition, Morris’ LSA formulation used elliptical coordinates to allow a bespoke solution scheme to be implemented, however, this limits the base flow (flow profile of the nozzle) to strictly elliptical geometries. However, elliptical jet base flow profiles are known to deviate significantly from a perfect elliptical profile downstream by switching axis and even bifurcating (Hussain & Husain 1989; Zeng et al. 2023). Shock-cell structures in elliptical jets also possess complex characteristics that can affect the development of the different flow structures (Menon & Skews 2010; Li et al. 2017; Nagata et al. 2022). Therefore, a model that can handle arbitrary geometry flow profiles is desirable in studying these jets. Furthermore, the existing elliptical jet LSA literature primarily considers how independent variables at a single jet operating condition affect the growth rates. For instance, Morris studied the effect of temperature ratio and displacement thickness on growth rate, only at a single Mach number and aspect ratio (Morris & Bhat 1995). However, different operating conditions may alter how sensitive growth rates are to a given independent variable. Therefore, the parameter space can be further explored by considering how characteristics of instability waves are affected by changing multiple parameters simultaneously, such as aspect ratio, Strouhal number, Mach number and shear layer thickness.

The existing literature on LSA of elliptical jets is primarily concerned with the most unstable coherent structures. However, neutrally stable waves are of significance in the study of some jet noise phenomena, particularly in jet screech. Screech is a jet noise phenomenon in which high-amplitude tones at discrete frequencies are produced over certain ranges of nozzle pressure ratios ($NPR$), caused by an aeroacoustic feedback process (Edgington-Mitchell 2019). The feedback process involves the downstream propagating Kelvin–Helmholtz (KH) instability wave (linearly unstable waves predicted by LSA) and has been shown to be closed by an upstream travelling neutral wave called the guided-jet mode (GJM) (Edgington-Mitchell et al. 2018; Gojon, Bogey & Mihaescu 2018; Nogueira et al. 2022b). The fundamental mechanism of screech is the same in both round and elliptical jets (Edgington-Mitchell et al. 2022), but the manifestation of mode staging is somewhat different in elliptical jets (Edgington-Mitchell, Honnery & Soria 2015). The differences may in part be due to differences in characteristics of the guided-jet and KH waves. However, the characteristics of the GJM for elliptical jets have not been explored in the literature to date. It is well established that the GJM only exists for a narrow range of frequencies, and these frequencies are dependent on the Mach number, temperature ratio, velocity profile and a number of other parameters (Tam & Hu 1989). One less well understood property of the GJM is its spatial dependence. The first results of Edgington-Mitchell et al. (2018) used a filtering in wavenumber space to reconstruct approximate amplitude profiles for the GJM; the mode was seen to exhibit a peak at some distance downstream, and decay away in either direction. This behaviour is not consistent with the linear stability approximation of this wave, which is that it is neutrally stable in the spatial sense. This behaviour was confirmed for a wider range of cases in Edgington-Mitchell et al. (2021) and Li et al. (2023). While care must always be taken in interpreting reconstructions based on filtered data, there is clear evidence that this wave, though predicted to be neutrally stable by linear stability analysis, exhibits some kind of decay behaviour in the physical jet. While this has been shown in multiple prior works, a definite explanation is yet to be provided. In the long-range resonance model of Mancinelli et al. (2021), this decay in space was hypothesised to be related to the ‘energisation’ of the GJM through a reflection mechanism. In the short-range resonance of Nogueira et al. (2022b), there is no requirement for the upstream wave to be neutral, and resonance is triggered by a pocket of absolute instability in the spatially periodic flow. Thus, the resulting global mode may be composed of waves that are spatially decaying, with their spatial growth rates only affecting the local spatiotemporal growth rate.

This work aims to predict the existence and characteristics (wavenumber, growth rate, spatial support) of instability waves in elliptical jets, and to further explore how operating conditions affect the behaviour of the supported instability waves. To that end, an elliptical jet velocity database will be explored and coupled with Reynolds-averaged Navier–Stokes (RANS) simulations so as to predict the overall characteristics of the mean flow of underexpanded elliptical jets. These will then be used as input to LSA to shed light on the salient features of the resonance loop. The current work is divided up as follows. In § 2, the experiments, numerical simulations and wave-modelling approaches are described. After that, the analysis of the experimental and simulation results are exposed in §§ 3 and 4, where an analysis of the different waves supported by the flow is also carried out. In § 5, the apparent disagreement between the spatially neutral LSA prediction and spatially decaying experimental result will be resolved. The paper is closed in § 6 with a revision of the main conclusions.

2. Methodology

2.1. Experimental set-up

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

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. (2013). 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 1987). 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 2023) solver. Details about the implementation and some verification cases for similar supersonic flow scenarios are given in Macia et al. (2019) and Heyns & Oxtoby (2014). The AUSM$^+$Up scheme proposed by Liou (2006) 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. (1995).

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. (2019) in the study of corrugated jets, and Nogueira & Edgington-Mitchell (2021) 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 (1989)), 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. (2019), 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. (2023). 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. 2019).

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 2009; Nogueira & Edgington-Mitchell 2021)

(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 2000), whereas the radial mesh is constructed from Chebyshev polynomials, which are then used to compute radial derivatives (Trefethen 2000). A transformation from Bayliss & Turkel (1992) 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. (2016); 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. (2018), Gojon et al. (2018); Edgington-Mitchell et al. (2021) and Nogueira & Edgington-Mitchell (2021)).

3. The elliptical jet base flow and the waves supported by it

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 (2022)) 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. (2013).

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 2021). 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 (1995) 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 (1995) 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. 2015). 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. 2018; Nogueira & Edgington-Mitchell 2021). 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. (2011) 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 (2021). 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. 2011; Edgington-Mitchell et al. 2021)

(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 1999; Edgington-Mitchell et al. 2021). 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. (2022b).

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. 2018, 2021; Nogueira & Edgington-Mitchell 2021), 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. 2022). 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. (2021), 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. Predictions from LSA

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 (1995) 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 (1995) 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. 2022a), 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. (2022a) 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 (1984) and Morris (2010)). 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. 2022b). While not directly related to the resonance model used in the present work, the long-range resonance model of Mancinelli et al. (2021) 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 (2017). 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 (1989) (see also Bogey (2021)) 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 1989). The structures in figure 18 resemble the ones computed for circular (Nogueira et al. 2022b) and twin (Nogueira & Edgington-Mitchell 2021) 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).

5. Spatial dependence of the GJM

As discussed in § 1, one of the most striking features of the GJM is that it only exists as a neutral wave over a discrete frequency range. An as-of-yet unexplained additional characteristic is that in experimental data, the amplitude envelope of the GJM does not match the expected shape for a neutral wave. An explanation may lie in the dependence of the neutral range on the thickness of the shear layer; Nogueira et al. (2022a) have previously shown for an $m=0$ circular jet that the upper and lower bounds of permitted Strouhal numbers for the GJM decrease with downstream distance (although the range of existence is still approximately the same). This hypothesis will now be explored with the current dataset. The vortex-sheet model (Bogey 2021) suggests that $m=1$ GJMs have a smaller existence region than $m=0$ GJMs, and that higher Mach numbers have a decreasing existence region. Therefore, one could expect that the $se_1$ GJM (which is reminiscent of $m=1$ disturbances) in elliptical jets will have a narrower frequency region than observed by Nogueira et al. (2022a) for the $m=0$ mode in the circular jet. The frequency region of existence of the $se_1$ GJM predicted by the LSA is shown in figure 20 for each $NPR$, computed at a precision of $\Delta St=0.01$.

Figure 20. Frequency range of existence of the $se_1$ GJM for $NPR=2.6$ (a,c) and $NPR=3.4$ (b,d).

It is evident that, at any axial location, the frequency range of existence predicted by the LSA is larger for the $NPR=2.6$ jet than the $NPR=3.4$ jet, consistent with the Mach number trend for circular jets observed by Bogey (2021). Furthermore, the maximum observed range of permitted Strouhal numbers at any axial location is $\Delta St = 0.09$ in the $NPR = 2.6$ jet and $\Delta St = 0.03$ in the $NPR = 3.4$ jet, both occurring when the shear layer is thickest. The upper and lower bounds of the frequency range decrease with downstream direction, also consistent with Nogueira et al. (2022a). This is likely due to the shear-layer thickness increasing with downstream direction, in line with the findings of Mancinelli et al. (2021), who demonstrated that increasing the shear-layer thickness decreases both the upper and lower bound of the frequency existence range. In the $NPR=3.4$ jet, the range of existence seems to oscillate whilst decreasing, especially for the RANS-based results. This is also consistent with Nogueira et al. (2022a), who observed a similar periodicity in axisymmetric jets, which was assumed to be a shock-cell effect. This is not observed in the $NPR=2.6$ flow, potentially due to the shock structure decaying more quickly in the axial direction. In both cases, the range of existence of the GJM computed using the RANS mean flows is usually smaller than the PIV-based ones, possibly due to the slower development of the shear layer in the downstream direction for the RANS (see figure 7).

With the relationship between streamwise position and the frequencies at which the GJM is supported by the flow established, an explanation for the observed amplitude envelope can be proposed. If the GJM was neutrally stable over the entire flow, it would not decay as it travels towards the nozzle. Instead, as demonstrated in figure 14, the guided jet mode peaks at some distance downstream, with a decay in either direction. In figure 20, the location at which the GJM attains maximum amplitude, as determined from figure 14, is plotted as a vertical dashed line. The experimentally determined screech frequency is plotted as a horizontal dashed line; as the POD modes capture fluctuations only at the screech frequency, the intersection of these lines is the relevant point of comparison with the LSA calculations. For both cases considered here, there is only a relatively narrow streamwise range where the GJM is predicted to be neutrally stable at the screech frequency. In both cases, the maximum observed amplitude of the GJM falls clearly within this narrow range for the PIV-based results. Predictions are slightly off for the RANS-based results, which is again due to the thinner shear layers predicted by the numerical model. This result is made clearer in figure 21, which shows the dispersion relation associated with the GJM computed at the streamwise position where the GJM amplitude derived from the POD is maximum. For both $NPR$, the screech frequency is within the bands of existence of this mode.

Figure 21. Dispersion relation of the $se_1$ GJM for $NPR=2.6$ (a) and $NPR=3.4$ (b) obtained where the POD mode is maximum.

The narrow streamwise range of existence shown in figure 20 may then be proposed as an explanation for the shape of the amplitude envelope observed in figure 14, in the context of the weakest-link model (Tam & Tanna 1982): downstream of this region, the GJM becomes damped, and any waves generated farther downstream will decay exponentially as they move upstream towards this region. Around the peak position, both KH and GJMs interact with the periodic shock-cells to trigger resonance. Upstream of this region, the guided jet mode cannot be supported as a discrete mode at this frequency; since the mode cannot suddenly vanish, it may be expected to become an acoustic wave travelling upstream, which are the only upstream waves supported by the flow at the screech frequency. In the context of the absolute stability model (Nogueira et al. 2022b), a similar mechanism may be proposed: the peak observed in the POD may represent a pocket of absolute instability in the flow (Chomaz, Huerre & Redekopp 1991; Pier 2008); downstream from that region, the spatiotemporal growth rate of the resonance may be negative, resulting in damped modes. Farther upstream, there is no GJM to trigger resonance, and the POD mode is simply dominated by the propagation of the modes. The long-range resonance scenario (Powell 1953; Mancinelli et al. 2019, 2021) provides equivalent explanations, with the peak of the POD mode representing the dominant position where the KH mode interacts with the shock-cell, and the upstream region being dominated by the propagation of the upstream mode. Embedded in all these explanations is the fact that all these resonance models consider a locally parallel (or, at best, spatially periodic) approach to approximate the upstream wave. Obviously, mean-flow development effects are not included in these approaches, so they can only capture part of the dynamics. The behaviour of the GJM as it travels in the upstream direction has not been studied in detail in the literature to the best of our knowledge but, given the acoustic-like nature of this mode, it may be hypothesised that, once it is not supported in its discrete form, it would simply become an acoustic wave, matching early descriptions of the screech phenomenon (Powell 1953). This hypothesis is corroborated by the fact that, even though the GJM is not a free stream acoustic wave, it bears several resemblances to acoustic waves travelling inside of the jet, being a manifestation of an acoustic wave generated in the region that experiences total reflection by the shear layer (Nogueira et al. 2024), a phenomenon that explains many of its key characteristics. Furthermore, it should be noted that the only upstream-travelling waves supported by the flow in the upstream region are acoustic waves that belong to the continuous spectrum (Gloor, Obrist & Kleiser 2013); as such, it may be expected that upstream waves in that region may be represented by a superposition of these acoustic waves.

6. Conclusions

In this paper, the dominant features of screeching elliptical jets were analysed using LSA. A PIV dataset for an elliptical aspect ratio 2 converging nozzle at two different Mach numbers was explored, and some characteristics of the flow could be identified. A complex shock-cell pattern was observed in the major axis, which included a low-speed streak that arises as to match the velocities at the centreline. The lower Mach number case was shown to experience axis switching at around five equivalent diameters away from the nozzle, which was also evident from the computed equivalent aspect ratio of the flow as function of the streamwise coordinate $x$. The experimental data was fitted to an elliptical jet profile, which allowed for the definition of approximated cross-plane velocity fields, which were necessary for the linear stability computations. The RANS simulations were also run for the same geometry to confirm the soundness of the fitting methodology and to provide alternative mean flows that were used to test the sensitivity of the results to details of the mean flow. Results show that the elliptical-jet approximation is accurate enough in the upstream region, but it fails to capture the overall flow shape farther downstream. It was also observed that, while the RANS managed to capture some dominant features of the flow (such as the first shock-cell position and the decay of the shocks), it failed to capture the position of axis switching and the variation in spacing of the shocks, requiring further tuning of the turbulence models to reproduce the exact same trends of the experiments.

Proper orthogonal decomposition was used to identify the wave-like structures involved in the screech loop. The method led to the identification of a clear pair of waves in the minor axis associated with most of the kinetic energy of the flow; no clear wave dominance was observed in the major axis, suggesting that the jets were resonating at a flapping mode. A spatial Fourier analysis of the travelling wave built using the POD modes allowed for the identification of the upstream and downstream components of the resonance loop, where their contrasting spatial support was highlighted. The characteristics of these different waves were then analysed in light of the KH and GJMs computed using spatial linear stability. The favourable match between model and experiments for the transverse velocity shapes, the frequencies of existence of the modes and their respective wavenumber confirmed that the $se_1$ mode is responsible for the flapping disturbances in elliptical jets, which is the preferred resonant stage at higher Mach numbers. It also confirms that LSA is able to capture the dominant characteristics of the waves responsible for screech in elliptical jets, despite being only a linear model.

One of the most interesting results of the current analysis is the fact that the GJM is only supported in a small range of streamwise positions for the flows at hand. The position where the amplitude of that wave is maximum in the POD is within the region where this mode is predicted in the stability analysis, but it ceases to be predicted at the screech frequency as it travels upstream – the linear stability model only predicts pure acoustic waves travelling at an angle in the negative region of the wavenumber space for regions closer to the nozzle. As there is no discrete mode to trigger the absolute instability in that region (Nogueira et al. 2022b), resonance may not be sustained. This description is somewhat similar to the presence of pockets of absolute instability in wake flows (Chomaz et al. 1991); in this case, this pocket of instability can only be present where the GJM is present in the flow. However, as the signature of upstream-travelling waves is non-negligible close to the nozzle, it is further hypothesised that such signature is simply due to the propagation of this wave in the upstream direction, which may become more and more similar to a pure acoustic wave as it travels upstream, given the lack of requirement for it to be a discrete mode. Further understanding of this phenomenon may be attained by means of spatial marching methods (Towne & Colonius 2015; Zhu & Towne 2023), which is considered to be outside of the scope of this work. Nevertheless, the current results are in support to the recent description of screech involving the KH wavepacket, the GJM and the shock-cell structure, which is now extended to elliptical jets.