2.1 Governing equations and numerical algorithm

The Navier–Stokes equations for single-phase, constant-density, incompressible, Newtonian fluid motion in an inertial reference frame are

(2.1a,b ) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}_{i}}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}}\boldsymbol{u}_{i}\boldsymbol{u}_{j}=-{\displaystyle \frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}}+\unicode[STIX]{x1D708}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\boldsymbol{u}_{i}}{\unicode[STIX]{x2202}x_{j}\unicode[STIX]{x2202}x_{j}}},\quad {\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}_{i}}{\unicode[STIX]{x2202}x_{i}}}=0. & & \displaystyle\end{eqnarray}$$

Here, $t$ , $\boldsymbol{u}_{i}(x,y,z)$ , $p(x,y,z)$ and $\unicode[STIX]{x1D708}$ are the time, velocity vector, pressure and kinematic viscosity of the fluid, respectively. For constant fluid density, the density may be combined with the pressure term.

In this paper, an unstructured, finite-volume algorithm developed by Mahesh, Constantinescu & Moin (Reference Mahesh, Constantinescu and Moin2004) is used to solve the Navier–Stokes equations (2.1). The spatial discretization emphasizes the simultaneous conservation of discrete first-order quantities (i.e. momentum) in addition to second-order quantities, such as kinetic energy. In other words, $\sum \boldsymbol{u}_{i}\unicode[STIX]{x2202}(\boldsymbol{u}_{i}\boldsymbol{u}_{j})/\unicode[STIX]{x2202}x_{j}$ over all control volumes only has contributions from the boundary elements. In this method, Cartesian velocities, $u_{i}$ , and pressure, $p$ , are stored at the control volume centroid. Additionally, face-normal velocities, $v_{n}$ , are stored separately at the centroids of the faces. The algorithm has been validated and used to simulate a variety of complex flows, including a gas turbine combustor (Mahesh et al. Reference Mahesh, Constantinescu and Moin2004), free jet (Babu & Mahesh Reference Babu and Mahesh2004) and transverse jets (Muppidi & Mahesh Reference Muppidi and Mahesh2005, Reference Muppidi and Mahesh2007, Reference Muppidi and Mahesh2008; Sau & Mahesh Reference Sau and Mahesh2007, Reference Sau and Mahesh2008; Iyer & Mahesh Reference Iyer and Mahesh2016; Regan & Mahesh Reference Regan and Mahesh2017).

A fractional-step (sometimes called predictor-corrector) method is used to solve the governing equations (2.1). Time is advanced explicitly using the Adams–Bashforth second-order scheme for the predictor velocities, $u_{i}^{\ast }$ , through the momentum equation using two previous time steps. The predicted velocities are then interpolated using second-order symmetric averaging to obtain the predicted face-normal velocities. A Poisson equation for pressure is then derived by taking the divergence of the momentum equation and satisfying continuity. This is used in a corrector step to project the solution onto a divergence-free velocity field. The Poisson equation is solved using the Algebraic Multi-Grid (AMG) solver in the HYPRE library (Falgout & Yang Reference Falgout and Yang2002). After solving for $p$ , $\boldsymbol{u}_{i}$ and $v_{n}$ are corrected using the pressure gradient.

2.2 Linear stability analysis

Modal linear stability analysis is the study of the dynamic response of a base state (i.e. base flow) subject to external perturbations (see Theofilis (Reference Theofilis2011) for a review). In this paper, the incompressible Navier–Stokes equations (2.1) are linearized about a base state $\overline{\boldsymbol{u}}_{i}(x,y,z)$ and $\overline{p}(x,y,z)$ , which can vary arbitrarily in space. The flow field is decomposed into a base state subject to a small $O(\unicode[STIX]{x1D700})$ perturbation ( $\tilde{\boldsymbol{u}}_{i}$ ). Subtracting the base flow equations yields the LNS equations

(2.2a,b ) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\boldsymbol{u}}_{i}}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}}\tilde{\boldsymbol{u}}_{i}\overline{\boldsymbol{u}}_{j}+{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}}\overline{\boldsymbol{u}}_{i}\tilde{\boldsymbol{u}}_{j}=-{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{p}}{\unicode[STIX]{x2202}x_{i}}}+\unicode[STIX]{x1D708}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\tilde{\boldsymbol{u}}_{i}}{\unicode[STIX]{x2202}x_{j}\unicode[STIX]{x2202}x_{j}}},\quad {\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\boldsymbol{u}}_{i}}{\unicode[STIX]{x2202}x_{i}}}=0. & & \displaystyle\end{eqnarray}$$

Note that the same numerical techniques are used to solve the LNS equations (2.2) and the Navier–Stokes equations (2.1). A molecular viscosity is used to perform linear stability analysis since the same was used to obtain the base flow obtained by DNS.

The LNS equations may be rewritten as a system of linear equations, where $\unicode[STIX]{x1D63C}$ is the LNS operator and $\tilde{\boldsymbol{u}}_{i}$ is the divergence-free velocity perturbation field. Solutions to the linear system of equations are of the form

(2.3) $$\begin{eqnarray}\displaystyle \tilde{\boldsymbol{u}}(x,y,z,t)=\mathop{\sum }_{j=1}^{N_{\unicode[STIX]{x1D714}}}\hat{\boldsymbol{u}}^{j}(x,y,z)\text{e}^{\unicode[STIX]{x1D714}_{j}t}+\text{c.c.}, & & \displaystyle\end{eqnarray}$$

where $N_{\unicode[STIX]{x1D714}}$ is the number of eigenvalues. Also, $\unicode[STIX]{x1D714}_{j}$ and $\hat{\boldsymbol{u}}^{j}$ can be complex. This defines $Re(\unicode[STIX]{x1D714})$ as the growth/damping rate and $\text{Im}(\unicode[STIX]{x1D714})$ as the temporal frequency of the complex velocity coefficient ( $\hat{\boldsymbol{u}}_{i}$ ). The system of equations transforms into a linear eigenvalue problem,

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D734}\hat{\unicode[STIX]{x1D650}}=\unicode[STIX]{x1D63C}\hat{\unicode[STIX]{x1D650}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{j}=\text{diag}(\unicode[STIX]{x1D734})_{j}$ is the $j$ th eigenvalue and $\hat{\unicode[STIX]{x1D650}}=(\hat{\boldsymbol{u}}^{1},\hat{\boldsymbol{u}}^{2},\ldots ,\hat{\boldsymbol{u}}^{N_{\unicode[STIX]{x1D714}}})$ is the matrix of eigenvectors.

2.2.1 Solutions of the LNS equations

For linear stability analysis, the size of the eigenvalue problem (2.4) can be $O(10^{6}{-}10^{8})$ . This makes solving the eigenvalue problem using direct methods very computationally expensive, often prohibitively so. Instead, an extension of the Arnoldi iteration method (Arnoldi Reference Arnoldi1951) called the implicitly restarted Arnoldi method (IRAM) is used, which is a matrix-free method. The present work uses the IRAM implemented in the P_ARPACK library (Lehoucq, Sorensen & Yang Reference Lehoucq, Sorensen and Yang1997) to efficiently calculate the leading (i.e. most unstable) eigenvalues and their associated eigenmodes.

A temporal exponential transformation of the eigenvalue spectrum is performed. The eigenvalue problem (2.4) is integrated over some time, $\unicode[STIX]{x1D70F}$ , which yields the exponential of the eigenvalue problem (2.4) which can be written as

(2.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D72E}\hat{\unicode[STIX]{x1D650}}=\unicode[STIX]{x1D63D}\hat{\unicode[STIX]{x1D650}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{j}=\text{diag}(\unicode[STIX]{x1D72E})_{j}$ . The matrix exponential $\unicode[STIX]{x1D63D}=\text{e}^{\unicode[STIX]{x1D63C}\unicode[STIX]{x1D70F}}$ is a time integration operator, which represents a numerical simulation of the LNS equations (2.2) over time $\unicode[STIX]{x1D70F}$ . This method is therefore described as a time-stepper method. Note that the eigenvectors, $\hat{\unicode[STIX]{x1D650}}$ , are the same between the two eigenvalue problems (2.4) and (2.5). However, the eigenvalues of the original problem (2.4) must be recovered using the following relationship:

(2.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{j}={\displaystyle \frac{1}{\unicode[STIX]{x1D70F}}}\ln \unicode[STIX]{x1D70E}_{j}. & & \displaystyle\end{eqnarray}$$

Readers are referred to Regan & Mahesh (Reference Regan and Mahesh2017) for an extensive validation of the solver for a variety of flow problems including the JICF.

2.3 Adjoint sensitivity analysis

The adjoint of a linear operator can be defined using the generalized Green’s theorem (Morse & Feshbach Reference Morse and Feshbach1953). In this paper, the continuous adjoint to the LNS equations are defined in the same way as in Hill (Reference Hill1995) and Giannetti & Luchini (Reference Giannetti and Luchini2007) (and similar to Barkley, Blackburn & Sherwin (Reference Barkley, Blackburn and Sherwin2008)), using the generalized Lagrange identity (Ince Reference Ince1926). The adjoint equations are

(2.7a,b ) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\boldsymbol{u}}_{i}^{\dagger }}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}}\tilde{\boldsymbol{u}}_{i}^{\dagger }\overline{\boldsymbol{u}}_{j}-\tilde{\boldsymbol{u}}_{j}^{\dagger }{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{i}}}\overline{\boldsymbol{u}}_{j}=-{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{p}^{\dagger }}{\unicode[STIX]{x2202}x_{i}}}-\unicode[STIX]{x1D708}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\tilde{\boldsymbol{u}}_{i}^{\dagger }}{\unicode[STIX]{x2202}x_{j}\unicode[STIX]{x2202}x_{j}}},\quad {\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\boldsymbol{u}}_{i}^{\dagger }}{\unicode[STIX]{x2202}x_{i}}}=0. & & \displaystyle\end{eqnarray}$$

Note the opposite sign on the viscous term, which defines that the adjoint equations must be solved backwards in time. The adjoint equations can also be rewritten as a system of linear equations, where $\unicode[STIX]{x1D63C}^{\dagger }$ is the adjoint LNS operator and $\tilde{\boldsymbol{u}}_{i}^{\dagger }$ is the adjoint to the velocity perturbation field. Similar to the direct problem, we assume non-trivial solutions to (2.7) of the form

(2.8) $$\begin{eqnarray}\displaystyle \tilde{\boldsymbol{u}}^{\dagger }(x,y,z,t)=\mathop{\sum }_{i=1}^{N_{\unicode[STIX]{x1D714}}}\hat{\boldsymbol{u}}^{\dagger }(x,y,z)\text{e}^{-\unicode[STIX]{x1D714}_{j}t}+\text{c.c}. & & \displaystyle\end{eqnarray}$$

Note the negative sign in front of $\unicode[STIX]{x1D714}$ , which allows for the eigenvalues from linear stability and adjoint sensitivity to have growth rates that correspond to their time integration directions (i.e. adjoint $Re(\unicode[STIX]{x1D714})>0$ corresponds to growth backwards in time). The adjoint systems of linear equations can now be simplified to an eigenvalue problem (similar to (2.9))

(2.9) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D734}\hat{\unicode[STIX]{x1D650}}^{\dagger }=\unicode[STIX]{x1D63C}^{\dagger }\hat{\unicode[STIX]{x1D650}}^{\dagger }, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{j}=\text{diag}(\unicode[STIX]{x1D734})_{j}$ is the $j$ th eigenvalue (coincident with the eigenvalue from linear stability analysis) and $\hat{\unicode[STIX]{x1D650}}^{\dagger }=(\hat{\boldsymbol{u}}^{\dagger ,1},\hat{\boldsymbol{u}}^{\dagger ,2},\ldots ,\hat{\boldsymbol{u}}^{\dagger ,N_{\unicode[STIX]{x1D714}}})$ is the matrix of adjoint eigenvectors.

Hill (Reference Hill1995) explains how $\tilde{\boldsymbol{u}}_{i}^{\dagger }$ , the adjoint velocity perturbation field, highlights optimal points in the flow where the largest response to unsteady point forcing will occur in its associated direct eigenmode counterpart. In the present work, adjoint sensitivity stability analysis is used in conjunction with linear stability analysis to determine flow regions that are most sensitive to point momentum forcing.

The same time-stepper method is implemented to efficiently compute the leading adjoint eigenvalues and their associated eigenvectors. The solver has been validated for obtaining adjoint modes in appendix A.

2.3.1 Wavemaker

This paper also discusses the receptivity of the JICF to spatially localized feedback. Due to the non-normality of the eigenvalue problem associated with the JICF, adjoint solutions alone cannot describe the whole picture. Therefore, the product for each $j$ th pair of direct and adjoint global modes is computed as

(2.10) $$\begin{eqnarray}\displaystyle W_{j}(x,y,z)={\displaystyle \frac{\Vert \hat{\boldsymbol{u}}^{j}\Vert \Vert \hat{\boldsymbol{u}}^{\dagger ,j}\Vert }{\text{max}(\Vert \hat{\boldsymbol{u}}^{j}\Vert \Vert \hat{\boldsymbol{u}}^{\dagger ,j}\Vert )}}, & & \displaystyle\end{eqnarray}$$

which determines the region where the eigenvalues of $\boldsymbol{A}$ are most sensitive to localized feedback (Giannetti & Luchini Reference Giannetti and Luchini2007) – also called ‘wavemaker’ regions. Locations where $W\approx 1$ are sensitive to localized feedback. The value of $W$ may be interpreted as a quantification of a possible change in the eigenvalues as a result of applied perturbations in the given region (Ilak et al. Reference Ilak, Schlatter, Bagheri and Henningson2012). Additionally, Giannetti & Luchini (Reference Giannetti and Luchini2007) have shown that the eigenvalues from linear stability analysis and adjoint sensitivity analysis are sensitive to domain size changes when values of the wavemaker, $W_{j}$ , are substantially different from zero at locations close to the domain boundaries. In wavemaker results that follow, all of the isocontours are displayed with a value of $0.01$ , and are spatially located far from the edges of the domain.

2.4 Optimal perturbation analysis

Non-modal stability analysis, or optimal perturbation analysis, is widely used in the literature for a variety of problems, as reviewed by Schmid (Reference Schmid2007). The traditional concept of Lyapunov stability does not coincide with non-modal stability analysis. Furthermore, the shape of the eigenmodes from optimal perturbation analysis and traditional asymptotic-time stability analysis can vary significantly, as they describe stability from different perspectives.

The energy is normalized with the initial energy when describing the transient growth

(2.11) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{E(\unicode[STIX]{x1D70F})}{E_{0}}}={\displaystyle \frac{(\tilde{\boldsymbol{u}}(\unicode[STIX]{x1D70F}),\tilde{\boldsymbol{u}}(\unicode[STIX]{x1D70F}))}{(\tilde{\boldsymbol{u}}_{0},\tilde{\boldsymbol{u}}_{0})}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}$ is the time scale over which the transient growth is optimized. This value of $\unicode[STIX]{x1D70F}$ is often smaller than the $\unicode[STIX]{x1D70F}$ used in linear stability analysis and adjoint sensitivity stability analysis in the previous sections. The perturbation energy may also be expressed in terms of the evolution operators $\unicode[STIX]{x1D63D}$ and $\unicode[STIX]{x1D63D}^{\dagger }$ defined above:

(2.12) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{E(\unicode[STIX]{x1D70F})}{E_{0}}}={\displaystyle \frac{(\unicode[STIX]{x1D63D}\tilde{\boldsymbol{u}}_{0},\unicode[STIX]{x1D63D}\tilde{\boldsymbol{u}}_{0})}{(\tilde{\boldsymbol{u}}_{0},\tilde{\boldsymbol{u}}_{0})}}={\displaystyle \frac{(\tilde{\boldsymbol{u}}_{0},\unicode[STIX]{x1D63D}^{\dagger }\unicode[STIX]{x1D63D}\tilde{\boldsymbol{u}}_{0})}{(\tilde{\boldsymbol{u}}_{0},\tilde{\boldsymbol{u}}_{0})}}. & & \displaystyle\end{eqnarray}$$

We are interested in the initial perturbations, $\tilde{\boldsymbol{u}}_{0}$ , that result in the largest transient growth. Equation (2.12) reveals that $\tilde{\boldsymbol{u}}_{0}$ is determined by the eigenvalues and eigenmodes of the operator $\unicode[STIX]{x1D63D}^{\dagger }\unicode[STIX]{x1D63D}$ . The eigenvalue problem may be expressed as

(2.13) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D726}\boldsymbol{U}^{\ast }=\unicode[STIX]{x1D63D}^{\dagger }\unicode[STIX]{x1D63D}\boldsymbol{U}^{\ast }, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D706}_{j}=\text{diag}(\unicode[STIX]{x1D726})_{j}$ is the $j$ th eigenvalue (i.e. growth factor) and $\unicode[STIX]{x1D650}^{\ast }=(\boldsymbol{u}^{\ast ,1},\boldsymbol{u}^{\ast ,2},\ldots ,\boldsymbol{u}^{\ast ,N_{\unicode[STIX]{x1D714}}})$ is the matrix of perturbation eigenmodes. The leading eigenmode offers the largest transient growth for the specified value of $\unicode[STIX]{x1D70F}$ , but sub-optimal eigenmodes often provide valuable insight towards other flow mechanisms that generate energy growth. In the discussions to follow, sub-optimal eigenmodes (i.e. perturbations) have smaller growth rates, but can become significant to the overall dynamics if their growth rates are near leading growth rate.

The size of the eigenvalue problem associated with optimal perturbation analysis is also computationally expensive to solve. The same method, the IRAM, is used to solve for the leading growth factors and the associated optimal perturbations. However, for optimal perturbation analysis, the eigenvalue problem is already formulated as a time-stepper method. The right-hand side of (2.13) is analogous to integrating a velocity perturbation forward some time $\unicode[STIX]{x1D70F}$ through the LNS equations, then backwards for time $\unicode[STIX]{x1D70F}$ through the adjoint equations. Therefore, the IRAM solves for the leading eigenvalues and eigenmodes without any further manipulation. This approach has been validated in appendix B.

2.5 Base flow generation

Linear stability, adjoint sensitivity and optimal perturbation analyses require a base flow around which the governing equations are linearized. For increasingly complex and globally unstable flows, obtaining a steady-state solution becomes difficult and computationally expensive. Hence, other approaches are being followed to solve for base flow to study more interesting and complex problems. Selective frequency damping (Åkervik et al. Reference Åkervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006) is one such approach to obtain a steady-state solution where a forcing term which acts as a temporal low-pass filter is added to the right-hand side of the governing equations. Some knowledge of the lowest unstable frequency is required for choosing the filter width. In order to converge to a steady solution, the filter cut-off frequency must be lower than that of all of the flow instabilities. Although this method lends itself to easy implementation, the computational cost is governed by the range of time scales. Additionally, selective frequency damping fails to dampen instabilities that are non-oscillatory, as shown by Vyazmina (Reference Vyazmina2010).

Another option is to use a turbulent mean flow as base state. Perhaps the best known example where linear stability analysis about the turbulent mean flow succeeds over the steady-state solution is the oscillating wake of a circular cylinder (Barkley Reference Barkley2006). The solutions about both base states agree at the onset of instability, but the steady-state base flow fails to capture the observed vortex shedding frequency far away from the bifurcation point. Recent studies by Turton, Tuckerman & Barkley (Reference Turton, Tuckerman and Barkley2015) and Tammisola & Juniper (Reference Tammisola and Juniper2016) used turbulent mean flow as a base state to examine linear stability around a turbulent mean flow. Barkley (Reference Barkley2006) and Turton et al. (Reference Turton, Tuckerman and Barkley2015) showed that performing linear stability analysis on a turbulent mean flow as base state results in eigenvalues with small real part and non-zero imaginary part.

Using turbulent mean flow as base state for the present JICF requires further discussion. Since a turbulent mean flow is a solution to the Reynolds-averaged Navier–Stokes equations, a nonlinear Reynolds stress term is effectively added to the governing equations when the base flow equations are used. This translates into a mode-dependent Reynolds stress being present in the eigenvalue problem. A scale-separation argument, first introduced by Crighton & Gaster (Reference Crighton and Gaster1976), and more recently discussed in the review by Jordan & Colonius (Reference Jordan and Colonius2013), can be used to justify when the mode-dependent Reynolds stress term is negligible. Only for the modes of interest (typically low frequency and large scale) must the Reynolds stress term be shown to be unimportant. Regan & Mahesh (Reference Regan and Mahesh2017) demonstrated that $St$ computed using the linearization around turbulent mean flow matches the results of DNS of Iyer & Mahesh (Reference Iyer and Mahesh2016). This proves empirically that using a turbulent mean flow as a base state for stability and sensitivity analyses can provide meaningful physical insight into flow dynamics. Additionally, the maximal value of Reynolds shear stress is negligible in the near field for both the mean flow fields (see appendix C), providing further justification for using the turbulent mean flows as base states.

Mantič-Lugo, Arratia & Gallaire (Reference Mantič-Lugo, Arratia and Gallaire2014) proposed a self-consistent model to obtain a base state identical to turbulent mean flow without requiring full nonlinear DNS. Their model requires knowledge of the most unstable mode and the saturation amplitude is determined by requiring the mean flow to be neutrally stable. The model was shown to give excellent results for flow over a circular cylinder compared to DNS in the range $50<Re<110$ . Their model, however, may not be directly applicable to a more complex problem such as JICF where the dynamics is dominated by multiple frequencies and $Re$ is high.

4.1 Upstream shear layer

The upstream shear-layer linear stability eigenmodes for both case $R2$ (figure 7 a) and case $R4$ (figure 8 a) were discussed in detail by Regan & Mahesh (Reference Regan and Mahesh2017). The main difference between the direct modes for each case is that for case $R2$ the mode originates near the jet exit plane, whereas for $R4$ the mode is elevated.

Figure 7.

The $R2$ upstream shear-layer linear stability (a) and adjoint sensitivity analyses (b) eigenmodes along with their associated wavemaker (c,d). Symmetry plane contours show the vertical velocity of the base flow $\overline{v}$ .

Figure 8.

Similar to figure 7, but for the $R4$ upstream shear layer.

The adjoint eigenmodes (figures 7 b and 8 b) show that the direct modes are most sensitive to $y$ -direction momentum forcing along the upstream side of the jet nozzle, near the jet exit. Interestingly, for $R2$ the wavemaker region (figure 7 c,d) is localized along the upstream side of the nozzle. Conversely, $R4$ (figure 8 c,d) is most sensitive to localized feedback along the entire upstream shear layer.

The wavemaker results are consistent with the notion that the upstream shear-layer region transitions from absolute to convective instability as $R$ changes from $2$ to $4$ . For  $R2$ , the region most sensitive to localized feedback is dominated by the formation of the upstream shear layer, which is in direct contrast to case $R4$ , which is sensitive to localized feedback along the entire upstream shear layer. The tonal nature of case $R2$ is due to the fact that the location where the shear layer begins to roll up (and its frequency dominates) is in the same location where the ‘wavemaker’ is the strongest. Case $R4$ is not only weaker, but the wavemaker region extends along the upstream shear layer.

4.2 Asymmetries in the CVP

Smith & Mungal (Reference Smith and Mungal1998) observed in their JICF experiments that asymmetries may form in the time-averaged CVP for $R>10$ . Similar observations were made by Getsinger et al. (Reference Getsinger, Gevorkyan, Smith and Karagozian2014) in their experiments who concluded that an absolutely unstable JICF ( $R2$ ) is less likely to exhibit asymmetric mean profiles compared to the weaker, convectively unstable JICF ( $R4$ ). The reason behind this behaviour of the JICF is not fully understood; specifically, the reason behind why there is a preferential direction in certain configurations.

Figure 9.

Similar to figure 7, but for the $R2$ left-leaning asymmetric eigenmodes. However, the linear stability mode (a) is highlighted using isocontour levels of $\pm 0.0001$ to highlight the streaks in the CVP.

Figure 10.

Similar to figure 7, but for the $R2$ right-leaning asymmetric eigenmodes. However, the linear stability mode (a) is highlighted using isocontour levels of $\pm 0.0001$ to highlight the streaks in the CVP.

Figure 11.

Similar to figure 7, but for the $R4$ left-leaning asymmetric eigenmodes.

Figure 12.

Similar to figure 7, but for the $R4$ right-leaning asymmetric eigenmodes.

In the present work, we observe significant asymmetries in some eigenmodes. The direct (a) and adjoint (b) eigenmodes in figures 9 and 10 for case $R2$ , and figures 11 and 12 for case $R4$ are left-leaning and right-leaning, respectively. Their corresponding wavemaker regions are biased towards each side as well. The adjoint eigenmodes (b) are most sensitive to vertical point momentum forcing in a similar way to the upstream shear layers, but with biases to each side. The wavemakers (c,d) are located along the CVP, directly behind the collapse of the jet potential core. By animating the linear stability analysis eigenmodes (not shown) it is seen that the eigenmodes for both cases rotate with the CVP.

It is important to examine if there are asymmetries present in the turbulent mean flow that are causing asymmetric eigenmodes. For case $R2$ , the turbulent mean flow was extended to an average over $170\,000$ samples. A contour plot of a $yz$ -plane at $x=1.33D$ is shown in figure 13 for the original base flow (figure 13 a) and the extended base flow (figure 13 b). Some small differences are visible between the two mean flows near the wall and in the middle of the CVP. Linear stability analysis was performed using the extended base flow, and the results are shown in figure 14. The eigenvalues from both simulations agree very well in both growth rate and $St$ . This provides evidence that the asymmetric eigenmodes are not a construct of an asymmetric mean flow. The turbulent mean flow for case $R4$ is also shown in figure 15 using contours on the $yz$ -plane at $x=D$ . The mean flow for case $R4$ does not suffer from asymmetries, and therefore does not require an extended statistics calculation for comparison.

Figure 13.

The velocity magnitude of two turbulent mean flows is shown in the $yz$ -plane at $x=1.33D$ . The first base flow (a) used 54 000 turbulent flow field samples, whereas the second (b) used 170 000 samples.

Figure 14.

Comparison of asymmetric mode eigenvalues after additional statistic samples. The upstream shear-layer eigenvalue is also shown for both calculations to provide context.

Figure 15.

The velocity magnitude of turbulent mean flow is shown in the $yz$ -plane at $x=D$ .

Linear stability analysis results for case $R2$ originate much closer to the jet nozzle exit compared to case $R4$ . The adjoint modes provide valuable information regarding the sensitivity of these asymmetric instabilities to $y$ -direction point momentum forcing. Note the spatial and temporal length scales that characterize the regions where asymmetric instabilities are most sensitive. For instance, adjoint sensitivity analysis results for case $R2$ (figure 9 b) show much longer length scales in the circumferential direction just below the jet nozzle exit when compared to case $R4$ (figure 11 b). This knowledge, in conjunction with the frequency information, provides valuable information regarding the best location and frequencies to excite asymmetries.

Growth rates from the linear stability and adjoint sensitivity analyses are often discussed in terms of their relative strength. We can compute the relative strength of the asymmetric eigenmodes for each case by comparing them to the strength of their respective upstream shear-layer growth rates. The difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{R2}^{\ast }$ between the growth rates of asymmetric eigenmodes and the upstream shear-layer eigenmodes for case $R2$ is in the range $0.042\leqslant \unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{R2}^{\ast }\leqslant 0.047$ . However, for the $R4$ case, the difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{R4}^{\ast }$ is in the range $0.014\leqslant \unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{R4}^{\ast }\leqslant 0.009$ . Notice $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{R2}^{\ast }>\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{R4}^{\ast }$ over their entire ranges, suggesting that asymmetric modes and sensitivity to experimental asymmetries are more significant for $R4$ than $R2$ , consistent with experimental results (Smith & Mungal Reference Smith and Mungal1998; Getsinger et al. Reference Getsinger, Gevorkyan, Smith and Karagozian2014).

4.3 Downstream of the jet exit

Figure 16 shows one of the pairs of downstream linear stability and adjoint sensitivity eigenmodes that have lower frequencies and longer length scales as compared to those previously discussed for case $R2$ . Additionally, they persist far downstream along the wall. Adjoint sensitivity analysis (figure 16 b) reveals sensitivities in the jet nozzle near the exit, but also in the region where the incoming cross-flow wraps around the jet nozzle exit, hinting at an increased sensitivity to perturbations from within the cross-flow boundary layer. The wavemaker region (figure 16 c,d) reveals some minor asymmetries in the sensitivity to localized feedback, but is largely symmetric between each side of the CVP.

Figure 16.

Similar to figure 7, but for a representative pair of $R2$ downstream eigenmodes.

Figure 17.

Similar to figure 7, but for the $R4$ lowest frequency eigenmodes.

Looking at case $R4$ , the lowest frequency pair of eigenmodes is shown in figure 17. The linear stability analysis eigenmode (figure 17 a) has larger spatial length scales than all previous eigenmodes, and also branches downwards towards the wall. There is a bias towards the right-hand side of the symmetry plane. However, this likely only implies that there is another set of eigenmodes that are biased towards the left-hand side of the symmetry plane, but with a growth rate smaller than that of the resolved eigenvalues. Adjoint sensitivity analysis (figure 17 b) shows that the largest sensitivity is not around the jet nozzle exit, but is upstream of the jet nozzle exit within the cross-flow boundary layer. Fric & Roshko (Reference Fric and Roshko1994) performed experiments that studied the wake vortices of the JICF by seeding the incoming cross-flow boundary layer and discovered that wake vortices could be visualized, leading to the conclusion that fluid inside the incoming cross-flow boundary layer travelled downstream to form wake vortices. Adjoint sensitivity analysis results add to these experimental results by showing there is a connection between perturbations in the cross-flow boundary layer and downstream of the jet exit near the wall. The wavemaker highlights the asymmetry, which implies that it is likely that a mirrored low-frequency pair of eigenmodes also exists.

4.4 Downstream shear layer

First shown by Regan & Mahesh (Reference Regan and Mahesh2017), the eigenvalue with the highest growth rate for case $R4$ is associated with the downstream shear layer and is shown in figure 18(a). As $R$ approaches infinity, the JICF becomes a free jet, and the upstream and downstream sides become indistinguishable. The magnitude of shear in the downstream shear layer is typically higher for case $R4$ compared to case $R2$ (see figure 45 d in appendix C). This makes the downstream shear layer dynamically more important for case $R4$ compared to case $R2$ . The linear stability mode is elevated from the jet nozzle exit and is located along the downstream side of the jet. This instability is most sensitive at the formation of the downstream shear layer. This region would be difficult to actuate in a control application, since it would most likely be invasive to the flow field.

Figure 18.

Similar to figure 7, but for the $R4$ leading downstream shear-layer eigenmodes (figure 6 b). Note that (a) is generated using the $80$ -million-element grid.

The wavemaker is located where the downstream shear layer forms. By building upon the previous analysis, figure 18(c,d) shows the localization of the wavemaker, which has a strong resemblance to the absolutely unstable upstream shear layer of case $R2$ (figure 7). This is yet another reason to identify this region of the downstream shear layer as absolutely unstable. Furthermore, an extension to higher values of $R$ would suggest that a critical value $R_{crit}$ exists, at which point the downstream shear-layer region becomes convectively unstable.

5.1 Optimal perturbations for case $R2$

Figure 19 suggests some overall conclusions from optimal perturbation analysis of case $R2$ . For short times (i.e. $\unicode[STIX]{x1D70F}<1/St_{up}$ ), asymmetric perturbations dominate where the CVP forms, with sub-optimal perturbations growing along the downstream shear layer. For $\unicode[STIX]{x1D70F}=0.8$ , there are sub-optimal modes that symmetrically perturb both the downstream shear layer and the CVP. Not until the characteristic time scale (i.e. $\unicode[STIX]{x1D70F}\approx 1/St_{up}$ ) does perturbing the downstream shear layer become optimal. On this time scale, it becomes clear from figure 19 that asymmetric perturbations of the CVP are sub-optimal. This gives rise to other sub-optimal perturbations of the upstream shear layer that quickly become significant for larger time scales. When optimal perturbations are considered for long time horizons (i.e.  $\unicode[STIX]{x1D70F}>1/St_{up}$ ), the modes result in growth along the upstream shear layer. Furthermore, sub-optimal perturbations move from the downstream shear layer to complex perturbations with higher circumferential wavenumbers along the upstream shear layer for $3.2\leqslant \unicode[STIX]{x1D70F}\leqslant 4.9$ . Finally, the group of perturbations that have the lowest growth factors include a series of hybrid perturbations that grow along both the upstream and downstream shear layers, as well as a series of increasing circumferential wavenumbers. The three time horizons are discussed in detail below.

5.1.1 Short time horizon

For the short time horizon, optimal perturbations take advantage of mechanisms that act much faster than the characteristic time scale $1/St_{up}=1.54$ to produce energy. First and foremost are the asymmetric perturbations, which dominate energy growth over the short time scales $0.4\leqslant \unicode[STIX]{x1D70F}\leqslant 0.8$ . The states of these perturbations at $\unicode[STIX]{x1D70F}=0.8$ are shown in figure 21. These perturbations ride along the CVP in the base flow as they propagate in helical fashion downstream. Not only are there pairs of asymmetric perturbations, but a series of increasing circumferential wavenumbers characterize the next few sub-optimal eigenmodes for $\unicode[STIX]{x1D70F}=0.4$ (not shown). The asymmetric perturbations originate at the jet nozzle exit and propagate downstream on either side of the CVP.

Figure 21.

Case $R2$ , short time horizon, $\unicode[STIX]{x1D70F}=0.8$ , final state of the first (a) and the second (b) leading asymmetric perturbations.

The next sub-optimal perturbations for $\unicode[STIX]{x1D70F}=0.8$ grow along the downstream shear layer, which are shown in figure 22. The evolution of the downstream shear layer perturbations is characteristic of the downstream shear-layer roll-up observed in DNS (figure 4 a). These perturbations originate within the nozzle on the downstream side as shown in figure 22. For $\unicode[STIX]{x1D70F}=0.4$ , the sub-optimal downstream shear-layer perturbation generates ${\approx}68\,\%$ of the energy growth compared to the asymmetric perturbation. Comparatively, when $\unicode[STIX]{x1D70F}=0.8$ , the sub-optimal downstream shear-layer perturbation (figure 22) generates ${\approx}73\,\%$ of the energy growth compared to the optimal asymmetric perturbation (figure 21). This highlights the fact that as $\unicode[STIX]{x1D70F}$ increases the downstream perturbations are becoming efficient at producing energy.

Figure 22.

Case $R2$ , short time horizon, $\unicode[STIX]{x1D70F}=0.8$ , origination (a) and final state (b) of the sub-optimal downstream shear-layer perturbation.

The least effective sub-optimal perturbations that generate the lowest energy growth for $0.4\leqslant \unicode[STIX]{x1D70F}\leqslant 0.8$ include a series of hybrid and higher wavenumber versions of the previously shown perturbations. For $\unicode[STIX]{x1D70F}=0.4$ , the least efficient perturbations are the hybrid perturbations that generate energy from both the CVP and the downstream shear layer. The results are similar for $\unicode[STIX]{x1D70F}=0.8$ , which have a series of hybrid CVP and downstream shear-layer perturbations shown in figure 23 with increasing wavenumbers and decreasing growth factors.

Figure 23.

Case $R2$ , short time horizon, $\unicode[STIX]{x1D70F}=0.8$ , final state of the first (a) and the second (b) sub-optimal asymmetric perturbations.

5.1.2 Characteristic time horizon

The characteristic-time-scale optimal perturbations take advantage of processes that act on times of the order of the characteristic time scale $1/St_{up}=1.54$ to increase energy. As $\unicode[STIX]{x1D70F}$ is increased beyond the short time scale, there is a bifurcation. The highest energy growth optimal modes change from asymmetric perturbations (figure 25) to perturbations that grow along the downstream shear layer (figure 24). Figure 19 clearly shows the significant drop in energy growth for the asymmetric perturbations, which are now sub-optimal for $\unicode[STIX]{x1D70F}=1.6$ .

Figure 24.

Case $R2$ , characteristic time horizon, $\unicode[STIX]{x1D70F}=1.6$ , origination (a) and final state (b) of the optimal downstream shear-layer perturbation.

Figure 25.

Case $R2$ , characteristic time horizon, $\unicode[STIX]{x1D70F}=1.6$ , initial (a) and final (b) state of a sub-optimal hybrid asymmetric downstream shear-layer perturbation. Note that there is a second perturbation that is mirrored across the $z=0$ plane.

Sub-optimal perturbations of the upstream shear layer are found for the first time in the optimal analysis at the characteristic time scale. The upstream shear-layer sub-optimal modes are shown in figure 26, arranged with decreasing energy growth and increasing circumferential wavenumbers. It is significant that higher wavenumber counterparts to the upstream shear-layer mode are also found. This shows that the upstream shear-layer region is becoming a significant energy growth opportunity since higher wavenumber downstream shear-layer and asymmetric perturbations are no longer recovered in the analysis.

Figure 26.

Case $R2$ , characteristic time horizon, $\unicode[STIX]{x1D70F}=1.6$ , initial (a) and final (b) state of the sub-optimal upstream shear-layer perturbation, along with the initial (c) and final (d) state of a higher circumferential wavenumber counterpart with a smaller growth factor.

5.1.3 Long time horizon

The long-time-scale optimal perturbations take advantage of processes that act on larger times than the characteristic time scale $1/St_{up}=1.54$ . Another bifurcation in optimal perturbation analysis for case $R2$ happens over the long time horizon. In the range $1.6\leqslant \unicode[STIX]{x1D70F}\leqslant 3.2$ the optimal perturbations change from acting along the downstream shear layer (figure 27) to the upstream shear layer (figure 28). Also, there are several higher wavenumber counterparts to the upstream shear-layer optimal perturbations (not shown).

Figure 27.

Case $R2$ , long time horizon, $\unicode[STIX]{x1D70F}=3.2$ , origination (a) and final state (b) of the sub-optimal downstream shear-layer optimal perturbation.

Figure 28.

Case $R2$ , long time horizon, $\unicode[STIX]{x1D70F}=3.2$ , the initial (a) and final (b) state of the upstream shear-layer optimal perturbation, along with the initial (c) and final (d) state of a lower growth factor perturbation with a higher circumferential wavenumber.

The least effective perturbations generate energy along the upstream and downstream shear layers simultaneously. For the longest time horizon that was studied for case $R2$ , the optimal perturbations again act along the upstream shear layer to generate the most energy growth, which are visualized in figure 29. This is consistent with the results from the linear stability and adjoint sensitivity analyses in § 4, which examined the stability and sensitivity at the asymptotic limit of time. Note that the first few sub-optimal perturbations are also upstream shear-layer modes with higher circumferential wavenumbers, one of which is shown in figures 29(c) and 29(d) for its initial and final states, respectively.

Figure 29.

Case $R2$ , long time horizon, $\unicode[STIX]{x1D70F}=4.9$ , initial (a) and final (b) state of the upstream shear-layer optimal perturbation, along the initial (c) and final (d) state of a lower growth factor perturbation with a higher circumferential wavenumber.

Figure 30 shows one of the remaining sub-optimal perturbations that are a set of hybrid upstream and downstream shear-layer modes. For these modes, there are no longer any purely downstream shear-layer perturbations. Instead, perturbations act on both the upstream and downstream side of the jet nozzle exit (figure 30 a) to generate growth along the jet trajectory (figure 30 b).

Figure 30.

Case $R2$ , long time horizon, $\unicode[STIX]{x1D70F}=4.9$ , origin (a) and final state (d) of the sub-optimal hybrid shear-layer perturbation.

5.2 Optimal perturbations for case $R4$

Similar to case $R2$ , some overall conclusions from optimal perturbation analysis for case $R4$ are highlighted using figure 20. Short-time-horizon optimal perturbations are dominated by growth along the downstream shear layer for $0.4\leqslant \unicode[STIX]{x1D70F}\leqslant 0.8$ . However, there are sub-optimal perturbations that grow along the upstream shear layer, as well as hybrid perturbations for both shear layers when $\unicode[STIX]{x1D70F}=0.8$ . For the characteristic time scale, the downstream shear-layer growth again dominates. However, a group of sub-optimal hybrid perturbations also have significant growth. For $\unicode[STIX]{x1D70F}=1.6$ , there are no longer purely upstream shear-layer sub-optimal modes, only hybrids that grow along the upstream and downstream shear layers. Over the long time horizons, it is shown that all optimal perturbations are some form of hybrid mode. When $\unicode[STIX]{x1D70F}=3.1$ , there is often a bias towards either the upstream or downstream shear layers. However, for $\unicode[STIX]{x1D70F}=4.7$ , the hybrid modes have less significant biases, with a more evenly spread perturbation across the upstream and downstream shear layers.

5.2.1 Short time horizon

The short-time-scale optimal perturbations take advantage of processes that act on times shorter than the characteristic time scale (i.e. less than $1/St_{up}=1.28$ ) to increase energy. Unlike case $R2$ , there are no asymmetric perturbation modes. Instead, the short-time-scale ( $0.4\leqslant \unicode[STIX]{x1D70F}\leqslant 0.8$ ) optimal modes, and the first few sub-optimal modes, generate energy growth as they propagate along the downstream shear layer. The group of downstream shear-layer modes for $\unicode[STIX]{x1D70F}=0.4$ are shown in figure 31, and for $\unicode[STIX]{x1D70F}=0.8$ in figure 32. For this short time, the perturbation modes do not have much time to advect along the base flow, which results in the focus on the downstream shear layer to generate growth.

Figure 31.

Case $R4$ , short time horizon, $\unicode[STIX]{x1D70F}=0.4$ , origination (a) and final state (b) of the leading downstream shear-layer perturbation.

Figure 32.

Case $R4$ , short time horizon, $\unicode[STIX]{x1D70F}=0.8$ , origination (a) and final state (b) of the leading downstream shear-layer perturbation.

The next few sub-optimal perturbation modes make use of the upstream shear layer to grow on top of the base flow. For $\unicode[STIX]{x1D70F}=0.8$ , figure 33(a) shows a sub-optimal perturbation mode with no circumferential wavenumber that generates the highest growth compared to its higher wavenumber counterparts (figure 33 b).

Figure 33.

Case $R4$ , short time horizon, $\unicode[STIX]{x1D70F}=0.8$ , initial (a) and final (b) state of sub-optimal upstream shear-layer perturbation, along with the initial (c) and final (d) state of a lower growth factor perturbation with a circumferential wavenumber.

5.2.2 Characteristic time horizon

Optimal perturbations on the characteristic time horizon take advantage of processes that act on times of the order of the characteristic time scale, $1/St_{up}=1.28$ , to increase energy. The optimal perturbations are again a group of modes that grow along the downstream shear layer. Figure 34 shows the leading optimal perturbation for $\unicode[STIX]{x1D70F}=1.6$ , along with the first few sub-optimal modes which also grow predominately along the downstream shear layer, but with higher circumferential wavenumbers.

Figure 34.

Case $R4$ , characteristic time horizon, $\unicode[STIX]{x1D70F}=1.6$ , initial (a) and final (b) state of the downstream shear-layer optimal perturbation, as well as the initial (c) and final (d) state of a sub-optimal perturbation with a smaller growth factor and higher circumferential wavenumber.

The remaining sub-optimal perturbation modes are all hybrid perturbation modes that show growth along upstream and downstream shear layers simultaneously. The hybrid perturbation with the largest growth factors is presented in figure 35. Hybrid modes have varying biases towards the upstream and downstream shear layers, and thus have slightly different evolutions. Interestingly, the final states of the perturbations at time $\unicode[STIX]{x1D70F}$ are optimized so that the upstream and downstream perturbations propagate along the shear layers and meet up at the collapse of the potential core.

Figure 35.

Case $R4$ , characteristic time horizon, $\unicode[STIX]{x1D70F}=1.6$ , origination (a) and final state (b) of a sub-optimal hybrid shear-layer perturbation.

5.2.3 Long time horizon

The long-time-scale optimal perturbations take advantage of processes that act on times longer than the characteristic time scale $1/St_{up}=1.28$ . The optimal perturbation for $\unicode[STIX]{x1D70F}=3.1$ continues to take advantage of the downstream shear layer as the only path to generate energy growth. The initial perturbation (figure 36) is elevated from the jet nozzle exit and grows significantly around the downstream shear layer as it travels further downstream. The sub-optimal hybrid perturbation with the largest growth factor is shown in figure 37; the other sub-optimal perturbations for $\unicode[STIX]{x1D70F}=3.1$ (not shown) again have different biases towards the upstream and downstream shear layers.

Figure 36.

Case $R4$ , long time horizon, $\unicode[STIX]{x1D70F}=3.1$ , origination (a) and final state (b) of the downstream shear-layer optimal perturbation.

Figure 37.

Case $R4$ , long time horizon, $\unicode[STIX]{x1D70F}=3.1$ , sub-optimal hybrid shear-layer perturbation with the largest growth factor.

For the longest time horizon of $\unicode[STIX]{x1D70F}=4.7$ for case $R4$ , all of the perturbations are shown to grow in energy along both upstream and downstream shear layers (figure 38). The initial perturbations of the hybrid modes (leading mode shown in figure 38 a) show biases towards both the shear layers, but their evolutions (leading mode shown in figure 38 b) are qualitatively similar. Therefore, for long time horizons for case $R4$ , the best way to perturb the base flow is to simultaneously perturb the upstream and downstream shear layers. This is consistent with the results from linear stability and adjoint sensitivity analyses in figure 6(b), which highlights that the upstream and downstream shear-layer instability modes have similar growth rates.

Figure 38.

Case $R4$ , long time horizon, $\unicode[STIX]{x1D70F}=4.7$ , the initial perturbation (a) and the associated final state (b) are shown for sub-optimal hybrid shear-layer perturbation with the largest growth factor.