Shooting method for linear inviscid bi-global stability analysis of non-axisymmetric jets

Abstract

The shooting method is commonly used to solve the linear parallel-flow stability problem for axisymmetric jets, i.e., a flow having one inhomogeneous direction. The present extension to two inhomogeneous directions – i.e., a bi-global stability problem – is motivated by inviscid non-axisymmetric jets. The azimuthal direction is Fourier transformed to obtain a set of coupled one-dimensional shooting problems that are solved by two-way integration from both radial boundaries – centreline and far field. The overall problem is formulated as one of iterative root-finding to match the solutions from the two integrations. The approach is validated against results from the well-established matrix method that discretizes the domain to obtain a matrix eigenvalue problem. We demonstrate very good agreement in two jet problems – an offset dual-stream jet, and a jet exiting from a nozzle with chevrons. A disadvantage of the shooting method is its sensitivity to the initial guess of the solution; however, this becomes an advantage when the need arises to track an eigensolution in a sweep over a problem parameter – say with increasing offset in the dual-stream jet, or with downstream distance from the nozzle exit. We demonstrate the performance of the shooting method in such tracking tasks.

Keywords

Shooting method bi-global stability Kelvin-Helmholtz instability inviscid non-axisymmetric jets

Introduction

Although axisymmetric (round) jets constitute a benchmark flow for their azimuthal homogeneity, practical prerogatives dictate the prevalence of non-axisymmetric jets in engines. The azimuthal inhomogeneity of such jets may be preferred, either for promoting mixing to reduce noise radiation, or for redirecting the noise away from the bottom sector of the jet. The former is exemplified by jets exiting from nozzles with chevrons,¹ and by jets from round nozzles having additional micro-jets impinging at their lip.² An instance of the latter is a dual-stream jet where the two streams are not coaxial, but instead have an offset between them such that the secondary potential core is thickened in the bottom sector.³

The linear Kelvin-Helmholtz (K-H) instability mode of the time-averaged flow field of turbulent jets is a useful model of their noise sources.⁴ A quasi-parallel flow assumption is often valid as the jet displays a slow streamwise spread. In case of axisymmetric jets, the consequent spatial stability analysis reduces to an eigenvalue problem involving an ordinary differential equation (ODE) in the radial coordinate, separately for every pair of frequency and azimuthal Fourier mode of fluctuation. However, for non-axisymmetric jets the problem becomes one of bi-global stability⁵ involving partial differentials in both the radial and azimuthal coordinates. This means that, in the Fourier azimuthal domain, although the linear stability equations still involve ordinary derivatives in the radial coordinates only, the various azimuthal Fourier modes of the eigenfunction are coupled for a particular frequency of perturbation. The latter problem has been solved using the matrix method,^6–9 wherein the system of coupled ODEs is converted into a matrix eigenvalue problem by suitable discretization of the radial domain. To be sure, there are many other instances^10–13 where the matrix version of the bi-global stability problem is solved in the physical polar coordinates or on a re-mapped Cartesian grid. The benefit of the azimuthal Fourier domain formulation is the ready simplification to the nominal axisymmetric jet, and the easy identification of the instability modes of the non-axisymmetric jets as continuation from their axisymmetric jet counterparts.

In this paper, we adopt the alternative approach of shooting.^14,15 In this method, the differential eigenvalue problem with boundary conditions is posed as an equivalent initial-value problem. The eigenvalue is guessed to start with. In the one-way shooting method, the eigenfunction is integrated starting from one boundary and proceeding towards the other. The satisfaction of the boundary condition thereat is obtained in an iterative manner by improving the guess of the parameters of the problem that include the eigenvalue. Basically, trajectories are ‘shot’ from one boundary in progressively more correct directions until one is found that hits the target at the other boundary. In the two-way shooting approach, the integration is started separately from both boundaries and approach each other at an intermediate point. The matching of the two eigenfunction solutions at this point is again achieved in an iterative manner.

The shooting method is preferred over the matrix method whenever (a) a single eigensolution is desired, and (b) a good initial guess is available for it. The bi-global stability problem will be seen to have multiple unstable K-H modes as solutions, and shooting may be used conveniently for ‘tracking’ these modes individually as some relevant condition (like axial station, Strouhal number, offset between the two jet streams, etc.) is varied. Another benefit of the shooting method is its reduced memory requirement compared to the matrix approach.

The shooting method has been used in the bi-global stability analysis of non-axisymmetric jets by Koshigoe et al.,^16,17 Morris et al.^18,19 and Gudmundsson.⁶ Here we report on some augmentation to the procedure for improved numerical stability and accelerated convergence; a preliminary version of this work appeared in Sohoni and Sinha.²⁰ We also provide validation of our shooting algorithm against the matrix method solution for an offset dual-stream jet and a jet exiting from a nozzle with chevrons. Moreover, we demonstrate how the shooting approach readily finds use in tracking of an eigenmode through incremental changes in some problem parameters, e.g., successive increments of the offset between the two streams of a dual-stream jet. This is a scenario for which the sensitivity of the shooting method to initial conditions makes it particularly well suited. The eigensolution for one parameter value is provided as the initial condition for the problem involving an incrementally different parameter value, thereby allowing rapid convergence of this solution.

Jets analyzed for validation

To motivate the development of the stability theory subsequently, we start by describing the kinds of non-axisymmetric jets that will serve as test cases for validation.

The first is an offset dual-stream jet where the ratio of the secondary to primary nozzle exit diameters ( $D_{s} / D_{p}$ ) is 1.7, the exit Mach numbers of the primary and secondary streams are respectively 1.5 and 0.9, and the entire jet is at ambient temperature such that there are no mean density variations. Figure 1(a) shows the contours of the mean axial velocity $\bar{u}$ of the jet at $x = 2 D_{p}$ with the offset C between the two streams being $0.1 D_{p}$ . The jet parameters were motivated by Murakami and Papamoschou,³ who presented experimentally-measured $\bar{u}$ profiles in the plane of symmetry in the cases of concentric and fully eccentric jets. Those profiles were fitted with the following velocity function,²¹ which provides a model whereby the effect of arbitrary variations of the offset may be investigated

\begin{array}{l} \bar{u} (x, y, z) = u_{c} (x) {(1 - h (x)) {\bar{u}}_{1} (x, y, z) + h (x) {\bar{u}}_{2} (x, y, z)}, \\ {\bar{u}}_{k} (x, y, z) = {\begin{array}{l} 1, if σ_{k} (y, z) < R_{k} (x), \\ \exp (- \frac{{(σ_{k} (y, z) - R_{k} (x))}^{2}}{δ_{k}^{2} (x)}), otherwise, \end{array} k \in {1, 2}, \\ σ_{1} (y, z) : = \sqrt{y^{2} + z^{2}}, σ_{2} (y, z) : = \sqrt{y^{2} + {(z + C)}^{2}} . \end{array}

(1)

Basically, two truncated Gaussian functions simulating the two streams are superposed, and the data of Murakami and Papamoschou³ is used to fit the parameters $u_{c}$ , h, R₁, R₂, $δ_{1}$ and $δ_{2}$ at each axial station (see Figure 1(c)). The offset C remains a free parameter that can be varied from 0 to $(D_{s} / D_{p} - 1) / 2$ to assess the effect of any intermediate offset between the concentric and fully eccentric extremes studied by Murakami and Papamoschou.³ Singh et al.²¹ applied the matrix approach of stability analysis to such jets with various offsets at a range of axial stations and perturbation frequencies; here we use the jet to validate the shooting approach. For later reference, all length dimensions in this jet are normalized by D_p, and the temporal angular frequency $ω$ of the perturbation analyzed is reported in terms of the Strouhal number $S t : = ω D_{p} / (2 π U_{p})$ , where $U_{p}$ is the primary jet exit velocity.

The second jet analyzed is the one exiting at Mach 0.9 from the SMC001 6-chevron nozzle designed and tested at NASA Glenn Research Center;¹ it was operated with a Mach 0.01 co-flow. The measured mean flow field is smoothed with fitting functions described by Sinha et al..⁸ The stability analysis is performed on the mean axial velocity field at $x = 0.5 D$ shown in Figure 1(a), where D is the nominal nozzle exit diameter. The density field is again assumed to be uniform. In this case, length dimensions are normalized by D, and $S t : = ω D / (2 π U)$ with U as the jet exit velocity.

In both jets, we ignore the cross-stream velocity fields, as well as any possible density/temperature variations. This is appropriate here since we are only setting out to validate the proposed bi-global shooting method with its matrix counterpart. The physics of the stability of these jets have been assessed in depth elsewhere.^8,21

Inviscid linear bi-global stability theory for non-axisymmetric jets

We use cylindrical coordinates $(x, r, θ)$ even though the formulation is for non-axisymmetric jets, as the reference nominal case is invariably an axisymmetric jet. Let $q : = {(u, v, w, p, ρ)}^{T}$ be the flow variable vector field involving the axial, radial and azimuthal velocity components, and pressure and density, respectively. The usual Reynolds decomposition is performed on q to obtain the time-averaged mean $\bar{q}$ and the residual fluctuations $q'$ . Assuming the mean flow to be locally parallel in x and stationary in time t, the linearity of the governing equations allow the following normal mode ansatz for $q'$

q' (x, r, θ, t) = {\tilde{q}}_{ω} (r, θ) e^{ι (α x - ω t)} + c . c .

(2)

In spatial stability analysis, $\tilde{q}$ is the eigenfunction corresponding to a specified real angular frequency $ω$ and the complex wave number $α$ is to be determined. The latter’s real part $α_{r}$ signifies the actual wave number that is inversely proportional to the phase speed $c_{p} : = ω / α_{r}$ ; its imaginary part $α_{i}$ corresponds to the decay rate (i.e., the negative of growth rate).

An inviscid analysis is typically warranted since the K-H instability is essentially an inviscid phenomenon and the jets under study have high Reynolds number. Then, substituting the above ansatz in the linearized compressible Euler equation obtains the usual compressible Rayleigh equation for pressure fluctuations. For the non-axisymmetric jets under analysis here, it takes the following bi-global form^7,11

\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial {\tilde{p}}_{ω}}{\partial r}) + \frac{1}{r^{2}} \frac{\partial^{2} {\tilde{p}}_{ω}}{\partial θ^{2}} - \bar{f} \frac{\partial {\tilde{p}}_{ω}}{\partial r} - \bar{g} \frac{1}{r^{2}} \frac{\partial {\tilde{p}}_{ω}}{\partial θ} - \bar{h} {\tilde{p}}_{ω} = 0

(3)

Here, the space-varying coefficients $\bar{f}, \bar{g}$ and $\bar{h}$ are functions of the mean flow field $\bar{q}$ as well as $α$ and $ω$

\bar{f} = \frac{2 α}{α \bar{u} - ω} \frac{\partial \bar{u}}{\partial r} + \frac{1}{\bar{ρ}} \frac{\partial \bar{ρ}}{\partial r}, \bar{g} = \frac{2 α}{α \bar{u} - ω} \frac{\partial \bar{u}}{\partial θ} + \frac{1}{\bar{ρ}} \frac{\partial \bar{ρ}}{\partial θ}, \bar{h} = α^{2} - \bar{ρ} {(α \bar{u} - ω)}^{2}

(4)

The mean pressure field is uniform in the free jets considered; the mean density field is also uniform in the particular jets analyzed here, but the corresponding gradient terms are retained for generality.

Instead of solving this problem in the physical $(r, θ)$ domain, we solve it in the Fourier azimuthal and physical radial domain. This facilitates subsequent specialization to axisymmetric jets, wherein the Fourier azimuthal modes of the solution are decoupled. Irrespective of the specific problem, the mean flow $\bar{q}$ and the pressure normal mode $\tilde{p}$ can be expanded in their respective azimuthal Fourier modes owing to periodicity in $θ$

\begin{array}{l} \bar{q} (θ) = \sum_{n = - \infty}^{\infty} {\hat{\bar{q}}}_{n} e^{ι n θ}, {\hat{\bar{q}}}_{n} : = \frac{1}{2 π} \int_{- π}^{π} \bar{q} (θ) e^{- ι n θ} d θ; \\ {\tilde{p}}_{ω} (θ) = \sum_{n = - \infty}^{\infty} {\hat{p}}_{ω, n} e^{ι n θ}, {\hat{p}}_{ω, n} : = \frac{1}{2 π} \int_{- π}^{π} {\tilde{p}}_{ω} (θ) e^{- ι n θ} d θ \end{array}

where we have omitted the

r -

dependence for notational compactness. For reference, Figure 1(b) shows the non-trivial azimuthal Fourier modes of the mean axial velocity fields of the jets studied here. The mean flow functions

\bar{f}, \bar{g}

and

\bar{h}

are also transformed similarly. With this, the Rayleigh equation becomes

\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial {\hat{p}}_{ω, m}}{\partial r}) - \frac{m^{2}}{r^{2}} {\hat{p}}_{ω, m} - \sum_{n = - \infty}^{\infty} {{\hat{\bar{f}}}_{n} \frac{\partial}{\partial r} + \frac{ι (m - n)}{r^{2}} {\hat{\bar{g}}}_{n} + {\hat{\bar{h}}}_{n}} {\hat{p}}_{ω, m - n} = 0, \forall m

(5)

Thus, for any non-axisymmetric jet, the eigensolutions are coupled in their azimuthal Fourier modes. On the other hand, in an axisymmetric jet only the zeroth azimuthal Fourier mode of the mean flow is non-trivial, thereby decoupling all the azimuthal modes of the eigensolution.

In the radial far field (i.e., $r \to \infty$ ), the base flow is uniform so that ${\hat{\bar{f}}}_{n} = {\hat{\bar{g}}}_{n} = 0$ for all azimuthal modes n, and ${\hat{\bar{h}}}_{n} = {\hat{\bar{h}}}_{0} δ_{n, 0}$ . This decouples equation (5) to the following modified Bessel equation thereat,²²

\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial {\hat{p}}_{ω, m}}{\partial r}) - (\frac{m^{2}}{r^{2}} + α^{2} - {\bar{ρ}}_{\infty} {(α {\bar{u}}_{\infty} - ω)}^{2}) {\hat{p}}_{ω, m} = 0, \forall m

(6)

Note that the $\bar{u}$ reduces to the uniform co-flow velocity ( ${\bar{u}}_{\infty}$ ) in the far field; correspondingly, $\bar{ρ}$ reduces to ${\bar{ρ}}_{\infty}$ . Unless there is a co-flowing jet, ${\bar{u}}_{\infty}$ will vanish in a laboratory setting. On the other hand, by replacing ${\bar{u}}_{\infty}$ with the mean centreline velocity ( ${\bar{u}}_{0}$ ) and ${\bar{ρ}}_{\infty}$ with the corresponding ${\bar{ρ}}_{0}$ , a similar simplification also occurs near the centreline (i.e., $r \to 0$ ) in the following circumstances. If the axial station under consideration is within the potential core (as is the case in this work), then the base flow is uniform at the centreline. If the jet has discrete rotational symmetry (e.g., a jet exiting from a nozzle with chevrons), then this centreline uniformity continues past the end of the potential core too. In case of isothermal offset dual-stream jets that are considered here, once we are past the primary potential core, we can shift the origin of the polar coordinates to the local peak of the mean axial velocity field to continue with the uniform base flow assumption at the modified centreline. However, in the case of a non-isothermal offset dual-stream jet, it may be impossible to find a suitable ‘centreline’ past the close of the potential core where the base flow can be assumed to be uniform. Barring this last extreme case, in all other non-axisymmetric jet problems, equation (5) reduces to the modified Bessel equation at both radial boundaries, yielding the boundary conditions

\begin{array}{l} {{\hat{p}}_{ω, m} |}_{r^{c}} \sim I_{m} (\sqrt{α^{2} - {\bar{ρ}}_{0} {(α {\bar{u}}_{0} - ω)}^{2}} r^{c}) = : I_{m} (λ^{c} r^{c}) \\ {{\hat{p}}_{ω, m} |}_{r^{f}} \sim K_{m} (\sqrt{α^{2} - {\bar{ρ}}_{\infty} {(α {\bar{u}}_{\infty} - ω)}^{2}} r^{f}) = : K_{m} (λ^{f} r^{f}) \end{array}

(7)

The boundary conditions are enforced at a small non-zero radius r^c (to avoid the centreline singularity), and at a very large radius r^f. Moreover, $I_{m}$ and $K_{m}$ are respectively the modified Bessel functions of first and second kind for azimuthal mode m.

As demonstrated here, all the equations are decoupled in $ω$ . So, for notational convenience, we will omit $ω$ in the subsequent development wherever it is obvious from the context.

Specialization to base flows with rotational symmetry

Chevron nozzles usually have the chevrons distributed uniformly around the circumference.¹ Similarly, nozzles with secondary micro-jets also typically have these devices deployed uniformly in azimuth.² Thus, the mean flow field in such jets exhibit an L–fold rotational symmetry, where L is the number of chevrons, micro-jets, etc. (see Figure 1(a) for an example). In such cases, the mean flow field presents a corresponding sparsity in the Fourier azimuthal domain:⁷

\bar{q} (r, θ) = \sum_{j = - \infty}^{\infty} {\hat{\bar{q}}}_{L j} (r) e^{ι L j θ}

(8)

This sparsity pattern induces a similar sparsity in the coefficient functions $\hat{\bar{f}}, \hat{\bar{g}}$ and $\hat{\bar{h}}$ . As a consequence, the Rayleigh problem of equation (5), with $ω$ omitted for notational convenience, becomes

\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial {\hat{p}}_{m}}{\partial r}) - \frac{m^{2}}{r^{2}} {\hat{p}}_{m} - \sum_{j = - \infty}^{\infty} {{\hat{\bar{f}}}_{L j} \frac{\partial}{\partial r} + \frac{ι (m - L j)}{r^{2}} {\hat{\bar{g}}}_{L j} + {\hat{\bar{h}}}_{L j}} {\hat{p}}_{m - L j} = 0, \forall m

(9)

Equation (9) indicates that, in the Fourier azimuthal domain the mth pressure azimuthal mode is coupled with the sparse set ${m + L j}_{j = - \infty}^{\infty}$ . To distinguish between the different (decoupled) eigenproblems, we define the azimuthal order M as the central azimuthal mode of the above coupled set. For example, the azimuthal order M = 0 eigenproblem couples the set of azimuthal modes ${L j}_{j = - \infty}^{\infty}$ ; the azimuthal order M = 1 problem involves ${1 + L j}_{j = - \infty}^{\infty}$ , and so on. Evidently, the unique azimuthal orders to solve for are $M \in {- ⌊ (L - 1) / 2 ⌋, \dots, - 1, 0, 1, \dots ⌊ L / 2 ⌋}$ , where $⌊ \cdot ⌋$ is the floor function. More details of these eigensolutions’ properties have been described by Sinha et al..⁸

Note that the above formulation reduces to the general non-axisymmetric jet case if we set L = 1, whereby all azimuthal modes are seen to be (densely) coupled. In this case, the only unique azimuthal order to solve for is M = 0.

Practical mean flow fields can be represented by a finite set of modes, say ${{\hat{\bar{q}}}_{L j}}_{j = - N}^{N}$ in equation (8). Owing to their nonlinearity, the corresponding mean flow functions $\hat{\bar{f}}, \hat{\bar{g}}$ and $\hat{\bar{h}}$ , have higher azimuthal complexity, say $\overset{⌣}{N} \geq N$ . Analogously, the eigensolutions will also converge with a finite coupled set of azimuthal modes in equation (9), say ${{\hat{p}}_{M + L j}}_{j = - S}^{S}$ , with $S \geq \overset{⌣}{N}$ . It will be noted that, if L is even, then the Nyquist azimuthal order $M = L / 2$ may be desired. In this case, symmetry considerations dictate that the coupled set of pressure eigenmodes to solve for is actually ${{\hat{p}}_{L (j + 1 / 2)}}_{j = - S}^{S - 1}$ .

Subsequently, it will be useful to identify and categorize the solutions of the eigenproblem of a particular azimuthal order M by the dominant azimuthal mode in the eigenfunction. Let us denote this dominant azimuthal mode number by $\tilde{m}$ . Clearly, $\tilde{m} \in {M + L j}_{j = - S}^{S}$ .

Specialization to base flows with mirror symmetry

Often, the jet nozzle geometry has a plane of symmetry, as in offset multi-stream jets or nozzles with symmetric chevrons, such that the resulting mean flow field also has a corresponding symmetry (see Figure 1 for two examples). In such cases, choosing the plane of symmetry as the $θ = 0$ reference, the mean axial velocity and density fields can be written as cosine series, with purely real azimuthal Fourier modes. Then, the azimuthal Fourier modes of the coefficient functions have the following symmetry properties

{\hat{\bar{f}}}_{- m} = {\hat{\bar{f}}}_{m}, {\hat{\bar{g}}}_{- m} = - {\hat{\bar{g}}}_{m}, {\hat{\bar{h}}}_{- m} = {\hat{\bar{h}}}_{m}

(10)

To deduce the consequent symmetries of the eigenfunctions, we replace m by – m and j by – j in equation (9), and use the above relations to obtain

\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial {\hat{p}}_{- m}}{\partial r}) - \frac{m^{2}}{r^{2}} {\hat{p}}_{- m} - \sum_{j = - \infty}^{\infty} {{\hat{\bar{f}}}_{L j} \frac{\partial}{\partial r} + \frac{ι (m - L j)}{r^{2}} {\hat{\bar{g}}}_{L j} + {\hat{\bar{h}}}_{L j}} {\hat{p}}_{- (m - L j)} = 0, \forall m

Comparing with equation (9), we observe that there is a one-to-one correspondence of the coefficients in the equations governing the positive and negative azimuthal mode counterparts in the eigensolutions.

We conclude that if L > 1 (i.e., if the flow has non-trivial rotational symmetry) and M is neither the axisymmetric nor the Nyquist azimuthal order, then the – M eigensolution can be retrieved from the $+ M$ one as

\begin{array}{l} M \in {0, L / 2} : α_{- M} = α_{+ M}, {\hat{p}}_{- (M + L j)} = β {\hat{p}}_{M + L j} \\ \forall j \in {- S, \dots, - 1, 0, 1, \dots, S}, \forall β \in ℂ \end{array}

(11)

with β reflecting the arbitrariness of the overall amplitude and phase of the eigenfunction in the linear homogeneous problem. In the remaining cases, viz. if L = 1 or L > 1 but M = 0 or the Nyquist azimuthal order, the negative and positive azimuthal mode counterparts of pressure eigenfunction are coupled in the same solution. Then, the preceding relation between the negative and positive pressure modes still hold, but now circularly within the same set. That is, we have

{\hat{p}}_{m} = β {\hat{p}}_{- m} = β (β {\hat{p}}_{m})

, which mandates that

β = \pm 1

. Thus, in this case, an eigenfunction can be either positive mirror symmetric or negative mirror symmetric; let us denote them by

{\hat{p}}^{+}

and

{\hat{p}}^{-}

, respectively. That is

\begin{matrix} M \in {0, L / 2} : {\begin{array}{l} Either \hat{p} \equiv {\hat{p}}^{+} s . t . {\hat{p}}_{- Lj}^{+} = {\hat{p}}_{Lj}^{+} \\ Or \hat{p} \equiv {\hat{p}}^{-} s . t . {\hat{p}}_{- Lj}^{-} = - {\hat{p}}_{Lj}^{-} \end{array}, \\ with j \in {\begin{array}{l} {0, 1, \dots, S} for M = 0, \\ \frac{1}{2}, \frac{3}{2}, \dots, (S - \frac{1}{2}) for M = L / 2. \end{array} \end{matrix}

(12)

Note that ${\hat{p}}_{0}^{-} \equiv 0$ . In the physical azimuthal domain, the positive and negative mirror-symmetric eigenfunctions have the following properties: ${\tilde{p}}^{+} (r, - θ) = {\tilde{p}}^{+} (r, θ)$ and ${\tilde{p}}^{-} (r, - θ) = - {\tilde{p}}^{-} (r, θ)$ ; this clearly justifies the nomenclature. In particular, we observe that ${\tilde{p}}^{-} (r, θ = 0) \equiv 0$ . Since the symmetry properties are unchanged if $θ$ is shifted by π, we also have that ${\tilde{p}}^{-} (r, θ = π) \equiv 0$ . Thus, the mirror plane is a nodal plane for the negative symmetric eigenfunction. We note here that a different, and less intuitive, terminology is also extant in the literature for these modes – cosinusoidal and sinusoidal, respectively.^11,12 Finally, it is observed that for mirror symmetric base flows, the unique azimuthal orders to solve for are $M \in {0, 1, \dots, ⌊ L / 2 ⌋}$ .

Starting from equation (9), the set of coupled Rayleigh equations governing the positive and negative mirror-symmetric eigensolutions of the M = 0 azimuthal order are

\begin{matrix} \frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial {\hat{p}}_{m}^{+}}{\partial r}) - (\frac{m^{2}}{r^{2}} + {\hat{\bar{f}}}_{0} \frac{\partial}{\partial r} + {\hat{\bar{h}}}_{0}) {\hat{p}}_{m}^{+} - \sum_{j = 1}^{N} ({\hat{\bar{f}}}_{L j} \frac{\partial}{\partial r} + {\hat{\bar{h}}}_{L j}) ({\hat{p}}_{| m - L j |}^{+} + {\hat{p}}_{m + L j}^{+}) \\ - \frac{ι}{r^{2}} \sum_{j = 1}^{N} {\hat{\bar{g}}}_{L j} {(m - L j) {\hat{p}}_{| m - L j |}^{+} - (m + L j) {\hat{p}}_{m + L j}^{+}} = 0, m \in {0, L, 2 L, \dots, S L} \end{matrix}

(13a)

\begin{array}{l} \frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial {\hat{p}}_{m}^{-}}{\partial r}) - (\frac{m^{2}}{r^{2}} + {\hat{\bar{f}}}_{0} \frac{\partial}{\partial r} + {\hat{\bar{h}}}_{0}) {\hat{p}}_{m}^{-} - \sum_{j = 1}^{N} ({\hat{\bar{f}}}_{L j} \frac{\partial}{\partial r} + {\hat{\bar{h}}}_{L j}) (sgn (m - L j) {\hat{p}}_{| m - L j |}^{-} + {\hat{p}}_{m + L j}^{-}) \\ - \frac{ι}{r^{2}} \sum_{j = 1}^{N} {\hat{\bar{g}}}_{L j} {| m - L j | {\hat{p}}_{| m - L j |}^{-} - (m + L j) {\hat{p}}_{m + L j}^{-}} = 0, m \in {L, 2 L, \dots, S L} \end{array}

(13b)

The advantage of explicitly enforcing the symmetry is evident; one needs to solve for only about half the total number of coupled azimuthal modes (the non-negative ones) of the eigenfunction. This not only halves the problem size, but also enforces the symmetries of the eigensolution exactly. Recall that M = 0 is the only azimuthal order to solve for in flows without rotational symmetry (i.e., when L = 1), making the computational gains particularly striking.

Matrix method

The matrix approach for this problem has been established over the past few years in a series of publications.^8,9,21 Hence, we will treat the results from this method as the ‘truth’, and validate the shooting approach proposed here with respect to them. We briefly outline the matrix procedure here; the details can be found in the above references.

The normal mode ansatz for the fluctuations described above is applied to the set of five linearized governing equations. The resulting eigenvalue problem (coupled in the azimuthal Fourier domain) is discretized using fourth-order central differences on a radial grid that is clustered close to the primary nozzle’s lip-line. The pole condition of Mohseni and Colonius²³ is applied at the centreline singularity, and the characteristic boundary condition of Thompson²⁴ is implemented at the far-field boundary ( $r = r^{f}$ ). This results in a generalized matrix eigenvalue problem that is solved using the ARPACK library²⁵ in a parallelized fashion. The implementation allows for the viscous effects to be retained or turned off; inviscid computations are made for the present comparison.

The main parameters for this algorithm are (a) the azimuthal modal complexity of the mean axial velocity N, (b) that of the eigenfunction solution S, (c) the radial location of the far-field boundary r^f, and (d) the number of points in the radial grid N_r. The first three parameters are shared with the shooting method too, but it will be obvious subsequently that they have subtle differences in their implications.

Shooting method for the bi-global stability problem

The shooting method is commonly used for solving two-point boundary value problems arising in one-dimensional linear stability problems.^14,15 To the knowledge of the authors, the only reported applications in bi-global stability problems are the works of Koshigoe et al.,^16,17 Morris et al.^18,19 and Gudmundsson;⁶ our formulation hews closest to the last reference. Here, we solve the bi-global Rayleigh equation (see equation (9)) for both serrated and offset dual-stream jets. Note that we need to determine an eigenvalue $α$ and the corresponding pressure eigenfunction ${{\hat{p}}_{ω, M + L j}}_{j = - S}^{S}$ for a particular choice of frequency $ω$ and azimuthal order M. Here, the Nyquist azimuthal order is not considered, as it requires special notation. The shooting method is most useful when only a few eigensolutions are desired, and good guesses are available for the corresponding eigenvalues. This happens to be the case for jet instability problems, as most often the unstable K-H mode is the only one of interest.

In this work, we employ a two-way shooting method,¹⁴ extending the one-way shooting approach described by Gudmundsson.⁶ The idea of a shooting method is to convert a two-point boundary value problem into an iterative initial value problem. We start with a guess of the eigensolution (that will be clarified below) at both the radial boundaries (centreline and far field), and shoot (i.e., integrate) them towards each other using equation (9). At a certain intermediate radial point rⁱ, say, the two eigenfunctions are compared through a cost function that is designed to be zero in case of a match. The process is necessarily iterative as the initial guesses have to be improved successively with the aim of zeroing the cost function till convergence is achieved.

In the one-way shooting approach, the guessed eigenfunction satisfying the applicable condition at one boundary is integrated to the other boundary, where the imposed condition is evaluated. The higher-order azimuthal modes of the eigenfunction have very small magnitudes at either radial boundary, making the evaluation of the match numerically inaccurate. The two-way method bestows greater numerical stability to the computations since rⁱ is chosen to be close to the peak of the eigenfunction.

The shooting starts with a guess of the eigenvalue $α$ that is to be refined progressively, and the integration of the pressure eigenfunction is initiated from both boundaries, i.e. from $r = r^{c}$ and $r = r^{f}$ (see equation (7)). We will denote these two eigenfunction solutions with superscripts ${(\cdot)}^{c}$ and ${(\cdot)}^{f}$ to respectively refer to the ‘centre’ and ‘far’ solutions integrated outward from the centreline and inward from the far field. The boundary condition in equation (7) indicates that, given the guess of $α$ , one knows the eigenfunction ${\hat{p}}_{m}$ as well as its first derivative at both r^c and r^f for all $2 S + 1$ coupled azimuthal modes involved. However, this is misleading since the relative amplitudes of the various ${\hat{p}}_{m}$ ’s must also be known – together the pressure modes have to satisfy the governing equation (9) over the entire r–domain. Thus, we refine the definition of the initial (or boundary) conditions presented in equation (7), and write them compactly as

{\hat{p}}_{m}^{ζ} (r^{ζ}) = Π_{m}^{ζ} B_{m} (λ^{ζ} r^{ζ}), (ζ, B) \in {(c, I), (f, K)}, m \in {M + L j}_{j = - S}^{S}

(14)

where

Π_{m}^{c}

and

Π_{m}^{f}

are complex scalar amplitudes to be determined. As the system is homogeneous, the overall amplitude of the eigenfunction is arbitrary so that one of the above scalars is a free variable in each solution. Recalling from § 3 that

\tilde{m}

is the expected dominant azimuthal mode in the eigenfunction, one choice is to set

Π_{\tilde{m}}^{c} = Π_{\tilde{m}}^{f} = 1

; we discuss a better choice subsequently. In any case, once these two values are fixed, all the remaining coupled azimuthal modes of the solution are fixed relative to each other. Thus, the vector of unknown parameters to be guessed in shooting is

G : = {(α, Π_{M - L S}^{c}, \dots, Π_{\tilde{m} - L}^{c}, Π_{\tilde{m} + L}^{c}, \dots, Π_{M + L S}^{c}, Π_{M - L S}^{f}, \dots, Π_{\tilde{m} - L}^{f}, Π_{\tilde{m} + L}^{f}, \dots, Π_{M + L S}^{f})}^{T}

(15)

These parameters will be uniquely determined in the correct solution that we are iterating towards.

To initiate the integration of the second-order ODE that is equation (9), we not only calculate the pressure eigenmodes at either boundary from equation (14) using the guessed parameter vector, but also their respective radial derivatives

{\frac{\partial {\hat{p}}_{m}^{ζ}}{\partial r} |}_{r = r^{ζ}} = Π_{m}^{ζ} λ^{ζ} B'_{m} (λ^{ζ} r^{ζ}), (ζ, B) \in {(c, I), (f, K)}, m \in {M + L j}_{j = - S}^{S}

(16)

Here, $B'$ denotes the derivative of the Bessel function B with respect to its argument. The integration from the boundaries towards the intermediate point rⁱ employs the variable-step-size Runge-Kutta solver ode45 in MATLAB®. For this, the second-order ODE is converted into a set of two coupled first-order ODEs.

Since $Π_{\tilde{m}}^{c}$ and $Π_{\tilde{m}}^{f}$ are chosen arbitrarily, the two solutions approaching from the two directions will not match at rⁱ in general (see Figure 2). However, in the correct eigenfunction, the ratios of each ${\hat{p}}_{m}^{c}$ to the corresponding ${\hat{p}}_{m}^{f}$ at rⁱ must be the same for all azimuthal modes m; the same ratios should also be maintained between $\partial {\hat{p}}_{m}^{c} / \partial r$ and $\partial {\hat{p}}_{m}^{f} / \partial r$ at rⁱ. We choose the value of ${\hat{p}}_{\tilde{m}}^{c} / {\hat{p}}_{\tilde{m}}^{f}$ at rⁱ as the common value against which to compare the ratios of all other modes and their radial derivatives at rⁱ. Thus, the set of cost functions to be zeroed are

Λ_{m} : = {({\hat{p}}_{m}^{c} - \frac{{\hat{p}}_{\tilde{m}}^{c}}{{\hat{p}}_{\tilde{m}}^{f}} {\hat{p}}_{m}^{f}) |}_{r = r^{i}}, Λ'_{m} : = {(\frac{\partial {\hat{p}}_{m}^{c}}{\partial r} - \frac{{\hat{p}}_{\tilde{m}}^{c}}{{\hat{p}}_{\tilde{m}}^{f}} \frac{\partial {\hat{p}}_{m}^{f}}{\partial r}) |}_{r = r^{i}}, m \in {M + L j}_{j = - S}^{S}

(17)

Note that the cost function $Λ_{\tilde{m}}$ is trivial. Thus, the non-trivial vector of cost functions is

Γ : = {(Λ_{M - S L}, \dots, Λ_{\tilde{m} - L}, Λ_{\tilde{m} + L}, \dots, Λ_{M + S L}, Λ_{M - S L}^{'}, \dots, Λ_{\tilde{m} - L}^{'}, Λ_{\tilde{m}}^{'}, Λ_{\tilde{m} + L}^{'}, \dots, Λ_{M + S L}^{'})}^{T}

(18)

It will be observed that, for a certain choice of S, the cost function vector above and the parameter vector G in equation (15) have the same number of elements, viz. $4 S + 1$ . In a mirror-symmetric problem for M = 0, the negative azimuthal modes’ amplitude factors are omitted from G , as are their cost functions. In this case the sizes of both G and $Γ$ vectors are $2 S + 1$ .

The multi-dimensional Newton-Raphson method is used as the iterative algorithm to find the parameter vector G that zeroes the cost function vector $Γ$ . Specifically, for a trial G , we evaluate $Γ$ at the intermediate radial point rⁱ, as well as the Jacobian $J$ of $Γ$ with respect to G at this point. Then the Newton step for the parameter vector $δ G$ is given by the solution of the linear system of equations $J δ G = - Γ$ . The problem is well-posed as G and $Γ$ are vectors of the same size, so that the Jacobian $J$ is a square matrix.

Unlike Gudmundsson⁶ who evaluated the Jacobian numerically using multi-dimensional finite differences, we calculate it analytically as described in Appendix 1. This requires additional quantities to be integrated from the boundaries along with the ‘center’ and ‘far’ pressure eigensolutions (see Appendix 1). For a given choice of azimuthal complexity S, the resulting size of the vector to be integrated separately from each boundary becomes $4 (S + 1) (2 S + 1)$ . In a mirror-symmetric problem for the M = 0 azimuthal order eigensolution, this size reduces to $2 (S + 1) (S + 2)$ . Although the calculation with the analytical Jacobian is more complex, it has the advantage of rapid convergence.

The foregoing shooting formulation fails in case of stable eigensolutions, since the solution becomes singular at the critical layer.^14,15 The standard workaround is to locally distort the integration path into the complex domain,⁶ but this is not implemented as of now. Thus, we limit our solutions to the unstable part of the eigenspectrum.

An advantage of the shooting approach is that the analytical boundary condition (see equation (7)) can usually be applied at a smaller outer radius than the corresponding characteristic boundary condition in the matrix approach. In theory, both are applicable wherever the jet mean flow gradient becomes zero. However, in the matrix approach, the backward difference approximation of derivatives at the boundary incur significant errors if the magnitude of the eigenfunction is non-negligible. This issue is well exemplified in one of the validation cases described subsequently.

The shooting method presented above can be used to analyze the stability of any non-axisymmetric jet or wake in a locally-parallel setting. For example, one could analyze rectangular, triangular or elliptic jets with this approach. Of course, the Fourier azimuthal parametrization may be more or less efficient depending on the particular problem. We demonstrate the method with two types of jet in this paper.

Initial guess of unknown parameter vector

The most subtle aspect of shooting is the initial guess of the unknown parameter vector G . In shooting-based stability analysis of flows with a single inhomogeneous direction (like axisymmetric jets), $G \equiv α$ – a scalar that is relatively easy to guess. In bi-global stability analysis, G is a high-dimensional vector that complicates its initial guess.

After much trial and error, the following heuristic approach was found to work reliably. Let us assume that an initial guess of the eigenvalue $α_{0}$ is available, may be from physical reasoning or otherwise. The various coupled pressure azimuthal modes of an eigenfunction are expected to be of comparable magnitude near their peak; else they would not need to be included in the solution in the first place. The pressure eigenfunction, and hence all the coupled ${\hat{p}}_{m}$ ’s, may be expected to peak close to the radial point of highest mean velocity gradient. For serrated jets, this is the nominal lipline. For dual stream jets, this can be either the inner lipline or the outer lipline – we will see subsequently that different eigensolutions are associated with these two maxima in the mean velocity gradient. We propose to set this radial point as the intermediate radial station rⁱ where the ‘center’ and ‘far’ shooting integrations terminate and must be matched. This way, the evaluation of the cost function $Γ$ and its Jacobian $J$ suffers from the least numerical inaccuracy.

Finally, in the absence of further information at the initiation, the individual pressure azimuthal modes of the eigenfunction are assumed to resemble Bessel functions not only at the boundaries but all the way throughout the shooting up to rⁱ, and all of them are assumed to attain a value of unity thereat. Thus, an initial guess of the complex scalar amplitudes is

Π_{m, 0}^{ζ} = \frac{1}{B_{m} (λ_{0}^{ζ} r^{i})}, (ζ, B) \in {(c, I), (f, K)}, m \in {M + L j}_{j = - S}^{S}

In the above, $λ_{0}^{c}$ and $λ_{0}^{f}$ are respectively the values of λ^c and λ^f evaluated at $α = α_{0}$ . Note that, $Π_{\tilde{m}}^{c} \equiv Π_{\tilde{m}, 0}^{c}$ and $Π_{\tilde{m}}^{f} \equiv Π_{\tilde{m}, 0}^{f}$ since these factors are not modified in the iterative shooting procedure.

The above discussion pertains to what we term ‘cold start’ of the shooting, as sketched in Figure 3(a). Here, one is interested in directly obtaining the final shooting solution with multiple coupled azimuthal modes without any prior knowledge of their approximate initial values at the boundaries. This is inherently difficult as the shooting algorithm is extremely sensitive to the initial guess. This problem may be mitigated in ‘warm start’ sketched in Figure 3(b), where we first obtain a solution with a few coupled azimuthal modes (may be just with $m = \tilde{m}$ , the expected dominant azimuthal mode), and then progressively add more and more azimuthal modes, always using the previous solution to initiate the common azimuthal modes. Referring to Figure 3(b), suppose that the $\tilde{m} = 0$ solution of an L = 1, M = 0 mirror-symmetric eigenproblem is available for the case of $S = 2 = : S_{0}$ , and a more refined solution is desired with $S > S_{0}$ . Then, the available solution’s $Π_{1}^{c}$ and $Π_{2}^{c}$ not only form the initial guesses $Π_{1, 0}^{c}$ and $Π_{2, 0}^{c}$ , respectively, but the latter is also used to obtain the initial guesses $Π_{3, 0}^{c}, Π_{4, 0}^{c}, \dots$ per the following general expression

Π_{m, 0}^{ζ} = \frac{Π_{M \pm L S_{0}, 0}^{ζ} B_{M \pm L S_{0}} (λ_{0}^{ζ} r^{i})}{B_{m} (λ_{0}^{ζ} r^{i})}, m \in {M \pm L j}_{j = S_{0} + 1}^{S}, (ζ, B) \in {(c, I), (f, K)}

Basically, we continue assuming that the new azimuthal modes of pressure to be included resemble Bessel functions and that they reach the same peak value (possibly different from unity) at rⁱ as the closest highest-order azimuthal mode for which a reliable initial condition is known from the earlier solution (see Figure 3(b)).

Stepped initiation integration to avoid numerical issues

The various azimuthal modes of the pressure eigenfunction are uncoupled and become corresponding Bessel functions towards the radial boundaries of the integration domain. Higher-order azimuthal modes of the eigenfunction typically have very small values at the extremes of the radial domain, both towards the centreline as well as in the far field (see illustration in Figure 3, as well as later results). From these minuscule values, these azimuthal modes grow very rapidly towards the shear layer, where they contribute to the coupled solution. In fact, at the radial extremes these modes may take on values close to or even less than the integration tolerance of the variable-step-size Runge-Kutta solver, which is clearly inadmissible in the numerical solution. A way out is to initiate these higher-order modes from radial positions that are progressively closer to the shear layer. This stair-stepping of the initial radial position of the integration is termed ‘stepped-initiation’ here; it is illustrated in Figure 3.

We first identify a threshold below which the solution may be beset by numerical integration error. For instance, $10^{- 10}$ may be a useful choice if the integration tolerance is set to $10^{- 12}$ , say. As outlined in § 4, the initial guess value of the amplitude factor for an azimuthal mode of the pressure eigenfunction assumes that it behaves as the corresponding Bessel function throughout its radial span of integration. Therefore, if this scaled Bessel function falls below the chosen threshold value anywhere within the radial span, then we initiate the integration of the particular azimuthal mode from the threshold-crossing radial point. The stepped-initiation strategy applies to both the cold start and the warm start, as depicted in Figure 3. In particular, Figure 3(a) shows that the radial starting points of the third and fourth pressure azimuthal modes in the far field are $r_{3}^{f}$ and $r_{4}^{f}$ , respectively; both are less than the chosen far-field boundary radius r^f. Also, per Figure 3(b), the radial starting point of the fourth pressure azimuthal mode near the centreline is $r_{4}^{c} > r^{c}$ . Note that for successive warm starts, the radial starting points of the various azimuthal modes of the eigenfunction should be saved and reused.

Results of validation assays

We present results from the study undertaken to validate the shooting method against the matrix method. The code is made general enough to handle both the offset multi-stream jet (having a mirror symmetry) and the jet exiting from the chevron nozzle (having rotational symmetry in addition to mirror symmetry), as described in § 2.

Validation with offset dual-stream jet

The first validation case is the dual-stream jet with offset C = 0.1, whose mean flow at $x = 2$ is considered (see Figure 1(a)). In this case, L = 1 as the base flow in this problem does not posses any rotational symmetry. Thus, the only non-trivial azimuthal order to solve for is M = 0. Stability calculations are performed for $S t = 0.3$ . The azimuthal modal complexity of the mean axial velocity field is set to N = 3 for both shooting and matrix approaches, as this forms their common input (see Figure 1(b) for reference).

Before discussing the validation results, we describe the eigensolutions for this problem. Table 1 presents a part of the unstable eigenspectrum, in terms of the growth rate $- α_{i}$ and the phase speed c_p. The inner and outer shear layers of the dual-stream jet have respective inner and outer K-H instability modes associated with them.^9,21,26 Owing to the Mach 1.5 primary jet, the inner modes display supersonic phase speeds. Conversely, the outer modes have subsonic phase speeds, as the secondary jet is subsonic. Apparently, there are several unstable inner and outer eigenmodes (we only present a few of them here); their labelling scheme is clarified subsequently. At this frequency ( $S t = 0.3$ ), the outer modes are more unstable than the inner ones.

Table 1.

Offset jet: calculations with the matrix and shooting methods yield almost identical eigenvalues of the $\tilde{m} = 0, 1^{+}$ and 1^- inner and outer shear layer modes.

	Inner modes					Outer modes
		Matrix		Shooting			Matrix		Shooting
$\bar{m}$	S	$c_{p}$	$- α_{i}$	$c_{p}$	$- α_{i}$	S	$c_{p}$	$- α_{i}$	$c_{p}$	$- α_{i}$
0	6	1.33379	0.26432	1.33378	0.26434	19	0.61789	0.27704	0.61790	0.27705
1⁺	6	1.26020	0.29162	1.26020	0.29162	15	0.59902	0.40404	0.59902	0.40404
1⁻	6	1.26424	0.29440	1.26424	0.29440	15	0.60023	0.34762	0.60023	0.34762

Results are for $S t = 0.3$ perturbations at $x = 2$ in case of offset C = 0.1. The azimuthal complexity S of the converged solutions is also presented.

Figure 4 shows the real part of the pressure eigenfunctions corresponding to the six eigenvalues presented in Table 1. The ‘far’ shooting solution for an eigenfunction is adjusted with a complex scalar to match its ‘centre’ solution counterpart at the intermediate radial grid point. Subsequently, the eigenfunction is normalized to have an absolute maximum of unity over the $y - z$ domain; moreover, the phase of the complex function is set to vanish at the point of maximum. The normalization is consistently applied to the matrix method results also. The mean flow depicted in Figure 1(a) should be referred for orientation of the axes maintained consistent here. The inner modes can be seen to peak around the inner shear layer ( $r \approx 0.5$ ) and the outer modes are maximum around the outer shear layer ( $r \approx 0.85$ ), justifying their nomenclature. Recall that r is normalized by the primary nozzle exit diameter, and the diameter ratio of the two nozzles is 1.7.

The rationale of the $\tilde{m}$ –labelling, already introduced in § 3, is explained in the context of the inner mode eigenfunctions shown in Figure 4. The $\tilde{m} = 0$ eigenmode is the offset jet continuation of the m = 0 mode in the concentric jet; it is dominated by the axisymmetric mode, although other azimuthal modes also contribute. It has positive mirror symmetry, and thus may be labelled $0^{+}$ . However, this is unnecessary as there is no $\tilde{m} = 0^{-}$ eigenmode (recall that eigenmodes with negative mirror symmetry have a trivial m = 0 component). The $\tilde{m} = 1^{+}$ and $1^{-}$ modes are respectively the positive and negative mirror symmetric continuations of the $m = \pm 1$ modes in the concentric jet; they continue to be dominated by these clockwise and counter-clockwise first helical modes. The series continues with $\tilde{m} = 2^{+}, 2^{-}$ , and so on, but these are not presented in the interest of brevity. This labelling is greatly confounded in the case of the outer modes; one has to actually track them with gradually increasing offset for disambiguation, as pursued subsequently.

Table 1 demonstrates the numerical similarity of the eigenvalues calculated by the shooting and matrix methods; they are evidently identical up to the four significant digits. These results are converged with respect to the main convergence parameter – viz. the azimuthal modal complexity S of the eigenfunction. In the shooting method, another parameter is $\overset{⌣}{N}$ ( $\leq S$ ) – the azimuthal modal complexity of the mean flow functions $\hat{\bar{f}}, \hat{\bar{g}}$ and $\hat{\bar{h}}$ . To reduce the number of free parameters, we set $\overset{⌣}{N} = S$ in all calculations here. The other azimuthal complexity parameter – viz. that of the mean flow field itself – is held at N = 3 for all these calculations for the purpose of comparison, as it is the common input to the two methods. Table 1 also demonstrates the different azimuthal complexities of the inner and outer modes. The inner shear layer is almost axisymmetric, the centreline being set to the centre of the primary jet. The asymmetry due to the offset between the two streams mainly manifests in the outer shear layer, thereby necessitating more azimuthal modes for the representation of the corresponding eigensolutions.

The validation of the eigenfunctions is demonstrated in Figure 4; in fact, the contours from the two approaches are so similar in this representation that they are not overlaid. A further demonstration of the similarity of the results from the two methods appears in Figure 5, where we resort to the azimuthal Fourier domain and plot the coupled modes on a logarithmic scale as in Figure 2. The first few azimuthal modes of the pressure eigenfunction are shown for the $\tilde{m} = 0$ inner and outer modes. The matrix method results are overlaid on those from the shooting approach, and the differences are indeed very small. The inner mode having supersonic phase speed decays comparatively slowly towards the far field, which is associated with its greater efficiency of noise radiation.⁴

The strategy of stepped initiation of integration outlined in § 4 was crucial to the calculation of these eigensolutions. Figure 5 shows that the higher order azimuthal modes of the $\tilde{m} = 0$ outer mode pressure eigenfunction start out with very small values indeed at both the radial bounds of the integration domain. These were automatically initiated from successively inward radial points, wherever they approximately approached the threshold value of $10^{- 10}$ . In fact, none of the individual azimuthal modes of the solution are above this threshold at $r = 10$ , so that the inward integration was initiated from $r^{f} = 5$ for all the outer mode cases. Similar stepped initiation is also apparent in the inner mode towards the centreline. Since these solutions decay slowly towards the far boundary, we used $r^{f} = 10$ for the inner modes. Note that the boundary values of these pressure azimuthal modes are not precisely equal to the set threshold value, since that is adjusted during the shooting iterations.

The minor differences in the normalized eigenfunctions calculated by the two methods are further highlighted in Figure 6(a), which also extends the presentation to the $\tilde{m} = 1^{+}$ and $1^{-}$ eigenmodes. The outer mode eigenfunctions calculated by the two approaches differ by less than 1 part in 10,000; recall that the eigenfunctions are normalized to a maximum of unity in the (r, θ) domain. However, the inner modes display greater discrepancies, specifically towards the far boundary.

To determine the cause of the discrepancy in the inner mode eigenfunction results, the matrix calculation of these modes is repeated with a higher value of r^f, the ‘far-field’ radius where we apply characteristic boundary conditions.²⁴ The previous inner mode results have been obtained with r^f=10 in both the shooting and matrix methods. Without redoing the shooting calculations, we evaluate their discrepancy against the matrix results recalculated with r^f=20. Figure 6(b) demonstrates that this reduces the error drastically; the eigenvalues were found to remain unchanged. In the matrix approach, the central difference scheme implemented within the radial domain must change to one-sided finite difference at the far boundary. Apparently, the inner modes have sufficient amplitude at r=10 (see Figure 5), so that the one-sided difference thereat incurs significant errors. The shooting method fares much better in this regard as it implements an analytical boundary condition that only requires the mean flow to be uniform at the boundary.

Finally, we note the parameters of the calculations that were not discussed above. For the shooting method, the ‘centre’ solution calculation was started from $r^{c} = 0.01$ . The intermediate radial point rⁱ where the matching of the two solutions was evaluated was set to the inner lipline ( $r = 0.5$ ) for the inner shear layer eigenmode calculations, and to the outer lipline ( $r = 0.85$ ) for the outer mode calculations. The absolute tolerance for the ‘zeroing’ of the cost function was set as $10^{- 5}$ ; this should be compared against the order-unity peak of the dominant azimuthal mode of the eigenfunction. The absolute and relative tolerances of the Runge-Kutta solver ode45 of MATLAB^® were both set to $10^{- 12}$ ; its variable step-size algorithm automatically used between 5000 and 10,000 steps in the radial domain. For the matrix method, the number of radial grid points N_r was 2500 in all the calculations, even when r^f was increased as in Figure 6(b); changing this did not affect the results materially.

Validation with single-stream chevron jet

The second validation case is the jet exiting from the 6-chevron nozzle analyzed at $x = 0.5$ (see Figure 1(a)). In this case there is a 6-fold rotational symmetry in the base flow; i.e., L = 6 in equation (9). The azimuthal modal complexity of the mean flow is N = 2 at this axial station. That is, only the m = 0, 6 and 12 azimuthal modes of the mean axial velocity are relevant (see Figure 1(b)). We solve the M = 0 eigenproblem for $S t = 0.3$ .

As described in § 3 and also discussed by Lajús et al.,¹¹ the eigenmodes of the six-chevron jet that are retrieved from the M = 0 problem can be categorized as $\tilde{m} \in {0, 6^{+}, 6^{-}, \dots}$ since L = 6. The $\tilde{m} = 0$ K-H eigenmode can be thought of as the continuation of the m = 0 K-H mode of the nominal round jet as the chevron progressively penetrates into the shear layer. The $\tilde{m} = 6^{+}$ and $\tilde{m} = 6^{-}$ K-H modes are respectively the positive and negative mirror-symmetric continuations of the $m = \pm 6$ helical K-H modes of the round jet. At the $x = 0.5$ station with $S t = 0.3$ , there are more unstable K-H modes; however, only $\tilde{m} = 0, 6^{+}$ and $6^{-}$ modes are presented here in the interest of brevity. Distinguishing between the various positive mirror-symmetric eigenfunctions (i.e., $\tilde{m} = 0, 6^{+}, 12^{+}$ , etc.) is not trivial. As Lajús et al.¹¹ have done, one has to track them from the round jet by continuously varying a chevron-impingement parameter. Similar effort is needed in case of the set of negative mirror symmetric eigenfunctions (i.e., $\tilde{m} = 6^{-}, 12^{-}$ , etc.). Here, we use the results from their work to label the eigenfunctions. Table 2 presents the eigenvalues of these modes calculated using the matrix and shooting methods; they are matched up to three significant digits.

Table 2.

Chevron jet: calculations with the matrix and shooting methods yield almost identical eigenvalues of the $\tilde{m} = 0, 6^{+}$ and 6⁻ K-H modes.

	Matrix			Shooting
$\tilde{m}$	$S$	$c_{p}$	$- α_{i}$	$S$	$c_{p}$	$- α_{i}$
0	13	0.7518	1.2378	8	0.7517	1.2384
6⁺	13	0.7040	0.7524	8	0.7040	0.7526
6⁻	13	0.4776	1.6213	11	0.4775	1.6211

Results are for $S t = 0.3$ , M = 0 perturbations at $x = 0.5$ . The azimuthal complexity S of the converged solutions is also presented.

The top row of the Figure 7 shows the real part of the pressure eigenfunction contours corresponding to the above mentioned modes. A 6-fold rotational symmetry is the underlying common feature displayed by all these M = 0 eigenfunctions. Also the positive and the negative mirror symmetry is evident for $\tilde{m} = {0, 6^{+}}$ and $\tilde{m} = 6^{-}$ solutions respectively.

Figure 7 also demonstrates the qualitative match of the eigenfunctions between the two approaches. This is further clarified in the bottom row that quantifies the difference in the Fourier azimuthal domain. The errors are less than 1 part in 1000, thereby validating the shooting approach. These errors should be considered in the context of the prevailing normalization of the eigenfunctions mentioned in the § 5; specifically, they reach a maximum value of unity in the $(r, θ)$ domain.

The stepped initiation of shooting is important in these calculations. It operates at the centreline and the far-field boundary from about m = 18 onwards in the all three eigenmodes. Without this artifice, it was impossible to obtain converged solutions.

These chevron jet results demonstrate that the shooting method is able to converge to the different unstable eigensolutions, depending on the initial guess. The initiation used the incremental ‘warm start’ strategy described in § 4, wherein we added one or two coupled azimuthal modes at a time to reach convergence in S. Throughout this process, the eigensolutions remained in the vicinity of their respective final (desired) values without veering off. The matrix approach is of course free from this issue as all the instability modes can be retrieved in one calculation. That the shooting approach is also able to pursue this task demonstrates the robustness of the implementation.

Tracking an eigensolution through parameter sweeps

In general, the bi-global stability problem is characterized by multiple unstable eigenmodes. As an example, the single-stream chevron jet possesses multiple concurrent instabilities.¹¹ Also, the offset round jet under investigation has unstable inner and outer K-H modes with $\tilde{m} = 0, 1, 2$ etc., all of which are simultaneous solutions of the same eigenproblem. Figure 4 showed that some of these instabilities are difficult to label unambiguously, and therefore to track with changes or sweeps of some problem parameters. Potential sweeping parameters of a bi-global stability problem are: the axial station x of the jet being studied, the Strouhal number St of the perturbation, the offset C between the two streams of a dual stream jet, the penetration of chevrons of a serrated nozzle, etc. In this section, we describe an algorithm that automates the process of sweeping through one of these parameters at a time to track a particular eigenmode, all the while ensuring convergence in the shooting parameters; it in turn calls the shooting program repeatedly.

To track a single eigenmode amongst several, we necessarily use the solution obtained with a particular sweep parameter value as a ‘warm start’ initial condition for the next increment of the parameter. On the one hand, the sensitivity of the shooting approach to initial conditions makes ‘cold start’ difficult; however, for the same reason, it is particularly efficient in case of warm start. Thus, the problem of tracking an eigensolution across a sweep of a parameter, where we repeatedly use warm starts to solve incremental problems, is particularly suited to the strengths of shooting.

There are multiple levels of complexity in the shooting method, all of which have their own convergence parameters. At the most fundamental level is the convergence with respect to the choice of integration tolerance of the variable-step size Runge-Kutta integration method. At the next higher level is the tolerance used for zeroing the shooting cost function vector $Γ$ . Then there are the choices of the far-field boundary radius, as well as the threshold to be set for the stepped-initiation. All of these parameters are preset to tight/high enough values, as determined in earlier studies described in the previous section; they are not varied in the parameter sweeps being pursued in this section.

The convergence parameters (termed shooting parameters here) that are iterated automatically for convergence in the present parameter sweep are the azimuthal complexity of the mean flow N, the azimuthal complexity of the mean flow functions $\overset{⌣}{N}$ (i.e., of $\bar{f}, \bar{g}$ and $\bar{h}$ ), and the azimuthal complexity of the eigenfunction S. These shooting parameters vary widely during any parameter sweep, and cannot be set to their highest possible values throughout the sweep since they directly determine the computational cost. The convergence in the three Fourier azimuthal mode parameters N, $\overset{⌣}{N}$ and S is the Fourier counterpart to the grid convergence test performed in any computation. The only difference is that there are three of these that appear in our approach. Selecting different values of these is analogous to having different physical azimuthal grids for representing the mean flow and the fluctuations, which is never done in calculations in the physical domain. Thus, for the purpose of the parameter sweeps discussed in this section, we chose all these three Fourier parameters to be the same, and here we will refer to them collectively as the azimuthal complexity S.

The S-convergence criterion is the closeness of the complex eigenvalue $α$ for two consecutive S values. That is, we choose the final S in a sequence of shooting evaluations such that $| {α |}_{S - 1} - {α |}_{S} | > ε_{α}$ and $| {α |}_{S} - {α |}_{S + 1} | \leq ε_{α}$ for a preset absolute tolerance $ε_{α}$ . The choice of this tolerance is dictated by two considerations: (1) the desired accuracy of the eigensolution, and (2) the proximity of the eigenmode being tracked to other simultaneous eigenmodes of the flow.

In general, the azimuthal complexity of the mean flow N increases with increasing offset C between the streams in a dual stream jet. It is expected that the azimuthal complexity of the eigenfunction S will also have the same trend with C. Hence, our sweeping program is set up to evaluate successively higher values of S for convergence in a sweep over increasing C. Although not pursued here, a sweep over increasing Strouhal numbers of perturbation will also necessitate evaluating increasing S values successively. On the other hand, the azimuthal complexity of the mean flow decreases with increasing axial position x downstream of the nozzle exit due to jet spread. For example, it is well known that jets issuing from serrated nozzles become approximately round by the end of the potential core.¹ To address this, our algorithm evaluates successively lower values of S when sweeping over increasing x.

The philosophy of warm start assumes that the eigensolution does not change significantly with the chosen increment of the sweep parameter. If the change is actually too drastic then the minimization of the shooting cost function may fail to converge when initiated from the previous solution. Our tracking algorithm automatically recovers from such a failure by reducing the step size of the parameter.

Tracking eigenmodes with increasing offset between streams of a dual-stream jet

Here we present results obtained when tracking eigenmodes in a dual-stream jet with offset C (normalized by the primary nozzle exit diameter D_p) being swept from 0 to 0.3. Recall that the ratio of secondary to primary nozzle exit diameters is 1.7, so that the offset can be at most 0.35, corresponding to a fully eccentric configuration. The calculations are performed with the mean flow at x=1, and the perturbations are at St=0.3. In the interest of brevity, the results presented here are limited to the inner and outer shear-layer eigenmodes for $\tilde{m} = 0$ and 1⁺. The absolute tolerance for eigenvalue convergence was set to ε_α=0.001; a looser tolerance prevented the correct tracking of eigenmodes in some extreme cases that will be described below. The azimuthal complexity automatically chosen by the sweeping program started from S = 0 (i.e., a single azimuthal mode) in the decoupled concentric case, and maximal values of 8 and 23 were reached at C = 0.3 in the tracking of inner and outer modes, respectively.

Figure 8 shows the variation of the growth rates and phase speeds of the various eigenmodes, with the offset being incremented in steps of 0.05. The monotonic trends observed here are consistent with correct tracking of the solutions. The similarity of the $\tilde{m} = 0$ and $\tilde{m} = 1^{+}$ outer eigenvalues at small offsets is what mandated the choice of the tight ε_α mentioned above.

The proper tracking of these eigenmodes is further validated by the consistent evolution of the corresponding eigenfunctions presented in Figure 9. Specifically, the real part of the pressure component is shown, and their mirror symmetry is invoked to present one half of the fields only. The complex eigenfunctions are scaled consistently so that their maximum values are unity, and their phase is set to 0° where this maximum is reached. Note that the orientation of these figures is such that the inner shear layer thickens in the bottom sector and thins out at the top with increase of offset.

The $\tilde{m} = 0$ inner eigenmode displays a gradual weakening in the top sector with increasing offset between the two streams. The $\tilde{m} = 1^{+}$ inner eigenmode also presents a similar trend, with the top lobe steadily splitting into two side lobes and the opposite phased bottom lobe shrinking in azimuthal extent. The $\tilde{m} = 0$ outer eigenmode, on the other hand, shows a clear preference for the top sector with increasing C, which is the opposite of the trend displayed by the $\tilde{m} = 1^{+}$ outer eigenmode. Looking at the two outer eigenfunctions obtained at C = 0.3, one would be hard pressed to label them as $\tilde{m} = 0$ or $1^{+}$ , since neither present any resemblance with the axisymmetric or first helical modes. However, the successful tracking of these modes from the C = 0 condition clarifies their identities unambiguously. This is all the more remarkable due to the closeness of their eigenvalues, as noted in Figure 8 as well as at the bottom of each panel in Figure 9. These findings were first presented by the authors in,⁹ and linked to the far-field sound in;²¹ those results were obtained using the matrix approach of § 3 and not the shooting method presented here.

The $C = 0^{+}$ notation deserves discussion. In the concentric jet (i.e., when C = 0), the azimuthal modes are decoupled. However, whenever the offset increases even infinitesimally, these positive and negative modes must combine into a positive and a negative mirror-symmetric solution, as shown in § 3. Thus, to initiate the tracking of the $\tilde{m} = 1^{+}$ positive mirror-symmetric mode here, one must combine the $m = + 1$ and m = – 1 solutions of the concentric jet appropriately, as would be found in case of $C \to 0^{+}$ . This exercise is not needed for the $\tilde{m} = 0$ eigenmode.

Tracking eigenmodes at successively downstream stations in offset dual-stream jets

A jet spreads as it flows downstream, so that the growth rate of a constant-frequency perturbation decreases monotonically. Therefore, such a perturbation grows, saturates and then decays, forming an axially-extended wavepacket. This is a very useful model for the large-scale coherent structures implicated in the loudest component of turbulent mixing noise in shear layers.⁴ In the quasi-parallel flow stability paradigm, a jet is assumed to be locally parallel. Thus, it is necessary to unambiguously track an unstable eigenmode found at a certain axial station by analyzing mean flow fields at successive downstream stations. This is facilitated by the shooting method, since its warm start approach is ideally suited to tracking eigensolutions across the variation of a sweeping parameter – the axial station in this case.

Axial tracking may be confounded if two or more eigenmodes approach each other at any axial station. It will be recalled from Figure 8 that for $S t = 0.3$ perturbations at $x = 1$ , the $\tilde{m} = 0$ and $\tilde{m} = 1^{+}$ inner shear layer eigenvalues are closest to each other when the offset in the dual stream jet is C = 0.2; for the outer modes, this occurs with C = 0.05. Figure 10 shows that shooting is successful at axial tracking of the four $S t = 0.3$ eigenmodes starting from $x = 1$ in these most problematic cases; the axial step size was $Δ x = 0.25$ . The tracking of each eigenmode was halted just upstream of the axial station where it stabilized. The inner shear layer modes remain unstable much farther downstream compared to the outer modes, reflecting the faster spread of the outer shear layer compared to the inner one (see Figure 1(c)).

The correctness of the axial tracking is further verified in the corresponding eigenfunctions presented in Figure 11. The distinct features of the $\tilde{m} = 0$ and $\tilde{m} = 1^{+}$ eigenmodes are maintained throughout the axial sweep. The inner eigenfunctions, that are tracked farther downstream, display a distinct broadening of the extent of the non-trivial structures in the eigenfunctions owing to the spread of the inner shear layer. The physical significance of these observations were explored at length by Singh et al..²¹

Conclusions

We describe a shooting approach for solving the bi-global inviscid parallel-flow linear stability problem for non-axisymmetric jets. Such a stability problem defined at a particular axial station of the jet may be formulated in polar coordinates or Cartesian coordinates. To connect with the limiting axisymmetric jet problem, we formulated the method in the Fourier azimuthal domain. In an axisymmetric stability analysis, the eigenproblem is fully decoupled in the Fourier azimuthal domain; on the other hand, non-axisymmetry of the base flow results in coupling of the Fourier modes of the eigenfunction. This increases the computational complexity of the problem.

There are a few formulations of the shooting approach to bi-global stability analysis reported in the literature.^6,16–19 We implemented several important improvements to the algorithm:

We use a two-way shooting approach instead of the more common one-way method, which suffers from greater numerical inaccuracy.

For faster and more reliable convergence, we calculate the necessary Jacobian of the cost function analytically instead of computations based on finite-differences.

We explicitly account for the mirror symmetry of the problem, thereby halving the problem size in many cases.

A novel ‘stepped-initiation’ strategy is introduced to allow inclusion of higher-order azimuthal modes of the solution that are too small in magnitude for accurate numerical handling near the boundaries of the radial domain.

Apart from detailing the formulation, we provide extensive evidence of the validity of our approach applied to two non-axisymmetric jets – a dual-stream jet with offset between the two streams, and a jet issuing from a serrated nozzle. Stability results obtained using the shooting method are compared against the prevalent matrix approach to the same problem, and excellent agreement is demonstrated in all the tests. It is also encouraging that the shooting approach is able to calculate various instability modes that are simultaneous solutions of the eigenproblem (as in the chevron jet near the nozzle exit) by suitable setting of the initial guesses.

The shooting method is especially useful in tracking eigensolutions in parametric sweeps. The solution obtained with one set of parameters may be conveniently used as the initial condition for solving with another parameter set that is not too different from the first, thereby hastening convergence of the iterative procedure. We demonstrate the correct tracking of concurrent eigenmodes with (a) increasing offset between the two streams of the dual-stream jet, and (b) increasing axial distance from its nozzle exit. The tracking algorithm introduced here automatically ensures convergence in the Fourier azimuthal complexity of the eigenfunction.

In this paper, we have not discussed the advantages of shooting method over the matrix approach with regards to computational efficiency, due to differences in implementation. The matrix method has been implemented in Fortran, and it is parallelized. Calculations are run on tens of processors, and results are obtained within a few minutes on a cluster of Intel Xeon multi-core processors, depending on the complexity of the calculation and the number of eigenvalues desired. On the other hand, the shooting method is implemented in serial mode in MATLAB^®, and comparable problems take tens of minutes on a MacBook Air laptop. Thus, although precise comparison would require commensurate codes and architecture, one can still draw conclusions on the comparable time-efficiency of the shooting method vis-à-vis the matrix approach. Similar computational benefits of the shooting approach have been brought out by others too.¹⁶

The shooting approach to the bi-global stability problem propounded here requires much less computer memory than the matrix method. In the latter approach, the high-dimensional matrix operators, although sparse, must be available in memory; a costly workaround is to evaluate the effect of the matrix operation on a candidate solution by re-calculating the matrix entries repeatedly. In the shooting method, only a one-dimensional slice of the eigenfunction is needed to be in memory at any given time – the values of all coupled azimuthal modes (and their various derivatives) at a certain radius. Given the current state-of-the-art in computer memory, this is admittedly not a significant advantage in bi-global stability analysis. However, future extension of the shooting approach to tri-global problems, if realized, may evince the benefit dramatically; the prevalent matrix approach to such problems challenges the limits of memory management.

The shooting method is applicable to a wide variety of analyses in physics. To our knowledge, all these implementations are for one-dimensional problems. The extension of the method demonstrated here to two dimensions for jets is readily applicable to all these physical problems, as necessary.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Nikhil Sohoni

Appendix 1. Analytically evaluating the Jacobian of the shooting cost function

We describe the analytical evaluation of the Jacobian $J$ of the cost function vector $Γ$ with respect to the parameter vector G . Equation (17) shows that $Γ$ has two kinds of entries, viz, $Λ_{m}$ and $Λ'_{m}$ , for various azimuthal modes m. It will be recalled that $m \in {M + L j}_{j = - S}^{S}$ in general, whereas $m \in {L j}_{j = 0}^{S}$ for the M = 0 azimuthal order of eigenfunction in mirror symmetric base flows. On the other hand, equation (15) shows that G consists of $α, Π_{n}^{c}$ and $Π_{n}^{f}$ for various azimuthal modes n. In general $n \in {M + L j}_{j = - S}^{S} ∖ \tilde{m}$ , but $n \in {L j}_{j = 0}^{S} ∖ \tilde{m}$ in the special case of the M = 0 solution in mirror symmetric base flows. The evaluation of the Jacobian thus requires the following derivatives at $r = r^{i}$ :

\begin{array}{l} \frac{\partial Λ_{m}}{\partial α} = \frac{\partial {\hat{p}}_{m}^{c}}{\partial α} - \frac{{\hat{p}}_{m}^{f}}{{\hat{p}}_{\tilde{m}}^{f}} \frac{\partial {\hat{p}}_{\tilde{m}}^{c}}{\partial α} - \frac{{\hat{p}}_{\tilde{m}}^{c}}{{\hat{p}}_{\tilde{m}}^{f}} (\frac{\partial {\hat{p}}_{m}^{f}}{\partial α} - \frac{{\hat{p}}_{m}^{f}}{{\hat{p}}_{\tilde{m}}^{f}} \frac{\partial {\hat{p}}_{\tilde{m}}^{f}}{\partial α}) \\ \frac{\partial Λ'_{m}}{\partial α} = \frac{\partial^{2} {\hat{p}}_{m}^{c}}{\partial α \partial r} - \frac{1}{{\hat{p}}_{\tilde{m}}^{f}} \frac{\partial {\hat{p}}_{\tilde{m}}^{c}}{\partial α} \frac{\partial {\hat{p}}_{m}^{f}}{\partial r} - \frac{{\hat{p}}_{\tilde{m}}^{c}}{{\hat{p}}_{\tilde{m}}^{f}} (\frac{\partial^{2} {\hat{p}}_{m}^{f}}{\partial α \partial r} - \frac{1}{{\hat{p}}_{\tilde{m}}^{f}} \frac{\partial {\hat{p}}_{\tilde{m}}^{f}}{\partial α} \frac{\partial {\hat{p}}_{m}^{f}}{\partial r}) \\ \frac{\partial Λ_{m}}{\partial Π_{n}^{c}} = \frac{\partial {\hat{p}}_{m}^{c}}{\partial Π_{n}^{c}} - \frac{{\hat{p}}_{m}^{f}}{{\hat{p}}_{\tilde{m}}^{f}} \frac{\partial {\hat{p}}_{\tilde{m}}^{c}}{\partial Π_{n}^{c}}, \frac{\partial Λ_{m}}{\partial Π_{n}^{f}} = \frac{{\hat{p}}_{\tilde{m}}^{c}}{{\hat{p}}_{\tilde{m}}^{f}} (\frac{{\hat{p}}_{m}^{f}}{{\hat{p}}_{\tilde{m}}^{f}} \frac{\partial {\hat{p}}_{\tilde{m}}^{f}}{\partial Π_{n}^{f}} - \frac{\partial {\hat{p}}_{m}^{f}}{\partial Π_{n}^{f}}) \\ \frac{\partial Λ'_{m}}{\partial Π_{n}^{c}} = \frac{\partial^{2} {\hat{p}}_{m}^{c}}{\partial Π_{n}^{c} \partial r} - \frac{1}{{\hat{p}}_{\tilde{m}}^{f}} \frac{\partial {\hat{p}}_{\tilde{m}}^{c}}{\partial Π_{n}^{c}} \frac{\partial {\hat{p}}_{m}^{f}}{\partial r}, \frac{\partial Λ'_{m}}{\partial Π_{n}^{f}} = \frac{{\hat{p}}_{\tilde{m}}^{c}}{{\hat{p}}_{\tilde{m}}^{f}} (\frac{1}{{\hat{p}}_{\tilde{m}}^{f}} \frac{\partial {\hat{p}}_{m}^{f}}{\partial r} \frac{\partial {\hat{p}}_{\tilde{m}}^{f}}{\partial Π_{n}^{f}} - \frac{\partial^{2} {\hat{p}}_{m}^{f}}{\partial Π_{n}^{f} \partial r}) \end{array}

The above expressions further suggest that we require the following set of quantities at the intermediate radial point rⁱ in order to evaluate the Jacobian (19)

{{\hat{p}}_{m}^{ζ}, \frac{\partial {\hat{p}}_{m}^{ζ}}{\partial r}, \frac{\partial {\hat{p}}_{m}^{ζ}}{\partial α}, \frac{\partial^{2} {\hat{p}}_{m}^{ζ}}{\partial α \partial r}, \frac{\partial {\hat{p}}_{m}^{ζ}}{\partial Π_{n}^{ζ}}, \frac{\partial^{2} {\hat{p}}_{m}^{ζ}}{\partial Π_{n}^{ζ} \partial r}}, ζ \in {c, f}

In the above, Einstein summing convention does not apply to repeated ζ indices. The first two quantities are found by default in the integration of the eigenfunction itself. To determine the four other quantities, we need to integrate them from the two boundaries also. Thus, the set of quantities in equation (19) actually constitute the augmented vector to be integrated. The necessary initial conditions at the boundaries for these additional terms are found as the corresponding derivatives of the boundary conditions in equations (14) and (16)

\begin{array}{l} {\frac{\partial {\hat{p}}_{m}^{ζ}}{\partial α} |}_{r = r^{ζ}} = Π_{m}^{ζ} r^{ζ} \frac{\partial λ^{ζ}}{\partial α} B'_{m} (λ^{ζ} r^{ζ}), {\frac{\partial {\hat{p}}_{m}^{ζ}}{\partial Π_{n}^{ζ}} |}_{r = r^{ζ}} = B_{m} (λ^{ζ} r^{ζ}) δ_{m, n}, \\ {\frac{\partial^{2} {\hat{p}}_{m}^{ζ}}{\partial α \partial r} |}_{r = r^{ζ}} = Π_{m}^{ζ} \frac{\partial λ^{ζ}}{\partial α} {B'_{m} (λ^{ζ} r^{ζ}) + r^{ζ} λ^{ζ} B''_{m} (λ^{ζ} r^{ζ})}, {\frac{\partial^{2} {\hat{p}}_{m}^{ζ}}{\partial Π_{n}^{ζ} \partial r} |}_{r = r^{ζ}} = λ^{ζ} B'_{m} (λ^{ζ} r^{ζ}) δ_{m, n} . \end{array}

As before, the pair $(ζ, B)$ takes values in $(c, I)$ or $(f, K)$ .

We differentiate equation (9) with respect to the relevant entries of G (viz. α and $Π_{n}^{ς}$ ) to obtain the second-order ODEs governing the last four quantities in the augmented integration vector of equation (19) (20)

\begin{array}{l} \frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial^{2} {\hat{p}}_{m}^{ζ}}{\partial r \partial α}) - \frac{m^{2}}{r^{2}} \frac{\partial {\hat{p}}_{m}^{ζ}}{\partial α} - \sum_{j = - \infty}^{\infty} {{\hat{\bar{f}}}_{L j} \frac{\partial}{\partial r} + \frac{ι (m - L j)}{r^{2}} {\hat{\bar{g}}}_{L j} + {\hat{\bar{h}}}_{L j}} \frac{\partial {\hat{p}}_{m - L j}^{ζ}}{\partial α} \\ - \sum_{j = - \infty}^{\infty} {\frac{\partial {\hat{\bar{f}}}_{L j}}{\partial α} \frac{\partial}{\partial r} + \frac{ι (m - L j)}{r^{2}} \frac{\partial {\hat{\bar{g}}}_{L j}}{\partial α} + \frac{\partial {\hat{\bar{h}}}_{L j}}{\partial α}} {\hat{p}}_{m - L j}^{ζ} = 0 \end{array}

(21)

\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial^{2} {\hat{p}}_{m}^{ζ}}{\partial r \partial Π_{n}^{ζ}}) - \frac{m^{2}}{r^{2}} \frac{\partial {\hat{p}}_{m}^{ζ}}{\partial Π_{n}^{ζ}} - \sum_{j = - \infty}^{\infty} {{\hat{\bar{f}}}_{L j} \frac{\partial}{\partial r} + \frac{ι (m - L j)}{r^{2}} {\hat{\bar{g}}}_{L j} + {\hat{\bar{h}}}_{L j}} \frac{\partial {\hat{p}}_{m - L j}^{ζ}}{\partial Π_{n}^{ζ}} = 0.

References

Bridges

Brown

CA.

Parametric testing of chevrons on single flow hot jets. In: 10th AIAA/CEAS aeroacoustics conference, AIAA Paper 2824, 2004.

Arakeri

Krothapalli

Siddavaram

, et al. On the use of microjets to suppress turbulence in a Mach 0.9 axisymmetric jet. J Fluid Mech 2003; 490: 75–98.

Murakami

Papamoschou

Mean flow development in dual-stream compressible jets. AIAA J 2002; 40: 1131–1137.

Jordan

Colonius

Wave packets and turbulent jet noise. Annu Rev Fluid Mech 2013; 45: 173–195.

Theofilis

Advances in global linear instability analysis of nonparallel and three-dimensional flows. Prog Aerosp Sci 2003; 39: 249–315.

Gudmundsson

Instability wave models of turbulent jets from round and serrated nozzles. PhD Thesis, California Institute of Technology, 2010.

Gudmundsson

Colonius

Spatial stability analysis of chevron jet profiles. In: 13th AIAA/CEAS aeroacoustics conference, AIAA Paper 3599, 2007.

Sinha

Gudmundsson

Xia

, et al. Parabolized stability analysis of jets from serrated nozzles. J Fluid Mech 2016; 789: 36–63.

Sohoni

Sinha

Modelling noise sources in offset two-stream jets using linear stability theory. In: 2018 AIAA/CEAS aeroacoustics conference, AIAA Paper 3941, 2018.

10.

Koshigoe

Tubis

Wave structures in jets of arbitrary shape. I. Linear inviscid spatial instability analysis. Phys Fluids 1986; 29: 3982.

11.

Lajús

Jr Sinha

Cavalieri

AVG

, et al. Spatial stability analysis of subsonic corrugated jets. J Fluid Mech 2019; 876: 766–791.

12.

Morris

Miller

DG.

Wavelike structures in elliptic jets. In: 22nd AIAA aerospace sciences meeting and exhibit, AIAA Paper 399, 1984.

13.

Uzun

Alvi

Colonius

, et al. Spatial stability analysis of subsonic jets modified for low-frequency noise reduction. AIAA J 2015; 53: 2335–2358.

14.

Drazin

Reid

WH.

Hydrodynamic stability. 2nd ed. Cambridge: Cambridge University Press, 2004.

15.

Schmid

Henningson

DS.

Stability and transition in shear flows. New York: Springer, 2001.

16.

Koshigoe

Tubis

Wave structures in jets of arbitrary shape. II. Application of a generalized shooting method to linear instability analysis. Phys Fluids 1987; 30: 1715.

17.

Koshigoe

Gutmark

Schadow

, et al. Wave structures in jets of arbitrary shape. III. Triangular jets. Phys Fluids 1988; 31: 1410–1419.

18.

Baty

Morris

PJ.

The instability of jets of arbitrary exit geometry. Int J Numer Meth Fluids 1995; 21: 763–780.

19.

Morris

Bhat

TRS.

The spatial stability of compressible elliptic jets. Phys Fluids 1995; 7: 185–194.

20.

Sohoni

Sinha

Shooting method for linear inviscid bi-global stability analysis of non-axisymmetric jets. In: 25th AIAA/CEAS aeroacoustics conference, AIAA Paper 2756, 2019.

21.

Singh

Sohoni

, and A. and Sinha Modelling noise from offset two-stream jets using linear stability theory and Kirchhoff surface method. In: 25th AIAA/CEAS aeroacoustics conference, AIAA Paper 2473, 2019.

22.

Michalke

Survey on jet instability theory. Prog Aerosp Sci 1984; 21: 159–199.

23.

Mohseni

Colonius

Numerical treatment of polar coordinate singularities. J Comput Phys 2000; 157: 787–795.

24.

Thompson

KW.

Time dependent boundary conditions for hyperbolic systems. J Comput Phys 1987; 68: 1–24.

25.

Lehoucq

Sorensen

Yang

ARPACK users’ guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods. SIAM, 1998.

26.

Gloor

Obrist

Kleiser

Linear stability and acoustic characteristics of compressible, viscous, subsonic coaxial jet flow. Phys Fluids 2013; 25: 084102.