2.1. Computational domain and general equations

The computational domain, represented in figure 2, models a coaxial flow configuration, which is composed of two inlet regions, an inner and outer pipe, both having a distance $D$ between walls and length $5D$, i.e. $z_{min} = -5 D$. The computational domain has an extension of $z_{max} = 50 D$ and $r_{max} = 25 D$. The distance between the pipes is equal to $L$, measured from the inner face of the outer tube to the face of the inner jet.

Figure 2. Computational domain of the configuration of two concentric jets used in StabFem.

The governing equations of the flow within the domain are the incompressible Navier–Stokes equations. These are written in cylindrical coordinates $(r,\theta,z)$, which are made dimensionless by considering $D$ as the reference length scale and $W_{o,max}$ as the reference velocity scale, which is the maximum velocity in the outer pipe at $z=z_{min}$

(2.1a)$$\begin{gather} \frac{\partial \boldsymbol{U}}{\partial t} + \boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{U} ={-}\boldsymbol{\nabla} P + \boldsymbol{\nabla} \boldsymbol{\cdot} \tau(\boldsymbol{U} ), \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U} = 0, \end{gather}$$

(2.1b)$$\begin{gather}\text{with}\ \tau(\boldsymbol{U}) = \frac{1}{{{Re}}} (\boldsymbol{\nabla} \boldsymbol{U} + \boldsymbol{\nabla} \boldsymbol{U}^T), \quad {{Re}} = \frac{W_{o,max} D}{\nu}. \end{gather}$$

Here, the dimensionless velocity vector $\boldsymbol{U} = (U,V,W)$ is composed of the radial, azimuthal and axial components, $P$ is the dimensionless reduced pressure, the dynamic viscosity $\nu$ and the viscous stress tensor $\tau (\boldsymbol{U})$.

The incompressible Navier–Stokes equations (2.1) are complemented with the following boundary conditions:

(2.2a,b)\begin{equation} \boldsymbol{U} = (0, 0, W_{i}) \text{ on } \varGamma_{in,i}\quad \text{and}\quad \boldsymbol{U} = (0, 0, W_{o}) \text{ on } \varGamma_{in,o}, \end{equation}

where

(2.3a,b)\begin{align} W_{i}=\delta_u \tanh{ \left(b_i (1 - 2r)\right)} \quad\text{and}\quad W_{o}=\tanh \left[ b_o \left( 1 - \left\lvert \frac{2r- (R_{outer,1}+R_{outer,2})}{D} \right\rvert \right) \right]. \end{align}

The parameter $\delta _u$ corresponds to the velocity ratio between the two jets, defined as $\delta _u=W_{i,max}/W_{o,max}$, the volumetric flow rates of the inner and outer jets are defined as $\dot {V}_{i} = 2 {\rm \pi}\int _0^{R_{inner}} r W_i \,\text {d}r$ and $\dot {V}_{o} = 2 {\rm \pi}\int _{R_{outer,1}}^{R_{outer,2}} r W_o \,\text {d}r$, respectively. The parameters $b_o$ and $b_i$ represent the boundary layer thickness within the nozzle, which is fixed equal to $5$ (as in Canton et al. Reference Canton, Auteri and Carini2017). With this choice of parameters, the volumetric flow rate of the inner jet is a function of the inner-to-outer velocity $\dot {V}_{i} = 3.73 \delta _u$, whereas the flow rate of the outer jet is a function of the duct wall length separating the two jets $\dot {V}_{o} = 5.41 L$. There is a weak influence of the boundary layer thickness on the stability properties of the jet, and it is related to the vortex-shedding regime developed upstream the separation wall (more details may be found in Talamelli & Gavarini (Reference Talamelli and Gavarini2006)). Finally, a no-slip boundary condition is set on $\varGamma _{wall}$ and a stress-free ($(({1}/{Re} )\tau (\boldsymbol{U})-P) \boldsymbol {\cdot } \boldsymbol{n} ={\boldsymbol {0}}$) boundary condition is set on $\varGamma _{top}$ and $\varGamma _{out}$, as shown in figure 2. In the following, the Navier–Stokes equations (2.1) and the associated boundary conditions will be written symbolically under the form

(2.4)\begin{equation} \boldsymbol{B} \frac{\partial \boldsymbol{Q}}{\partial t} = \boldsymbol{F}(\boldsymbol{Q}, \boldsymbol{\eta}) \equiv \boldsymbol{L} \boldsymbol{Q} + \boldsymbol{N}(\boldsymbol{Q},\boldsymbol{Q}) + \boldsymbol{G} (\boldsymbol{Q},\boldsymbol{\eta}), \end{equation}

with the flow state vector $\boldsymbol{Q} = [ \boldsymbol{U},{P} ]^{\rm T}$, $\boldsymbol {\eta } = [{{Re}}, \delta _u ]^{\rm T}$ and the entries of the matrix $\boldsymbol{B}$ arise from rearranging (2.1). Such a form of the governing equations takes into account a linear dependency on the state variable $\boldsymbol{Q}$ through $\boldsymbol{L}$. And a quadratic dependency on the parameters and the state variable through operators $\boldsymbol{G}({\cdot }, {\cdot })$ and $\boldsymbol{N}({\cdot }, {\cdot })$.

2.2. Asymptotic stability

2.2.1. Linear stability analysis

In this study, the authors attempt to characterise the stable asymptotic state from the spectral properties of the Navier–Stokes equations (2.1). First, let us consider the stability of an axisymmetric steady-state solution named $\boldsymbol{Q}_0$, which will be also referred to as the trivial steady state. For that purpose, let us evaluate a solution of (2.1) in the neighbourhood of the trivial steady state, i.e. a perturbed state, as follows:

(2.5)\begin{equation} \boldsymbol{Q}(\boldsymbol{x}, t) = \boldsymbol{Q}_0(\boldsymbol{x}, t) + \varepsilon \hat{\boldsymbol{q}}(r, z) \exp({-\mathrm{i} (\omega t - m \theta)}), \end{equation}

where $\varepsilon \ll 1$, $\hat {\boldsymbol{q}} = [ \hat {\boldsymbol{u}},\hat {{p}} ]^{\rm T}$ is the perturbed state, $\omega$ is the complex frequency and $m$ is the azimuthal wavenumber. The next step consists of the characterisation of the dynamics of small-amplitude perturbations around this base flow by expanding it over the basis of linear eigenmodes (2.5). If there is a pair $[\mathrm {i} \omega _\ell, \hat {\boldsymbol{q}}_{\ell } ]$ with $\text {Im} (\omega _\ell ) > 0$ (respectively the spectrum is contained in the half of the complex plane with negative real part) there exists a basin of attraction in the phase space where the trivial steady state $\boldsymbol{Q}_0$ is unstable (respectively stable) (Kapitula & Promislow Reference Kapitula and Promislow2013). The eigenpair $[\mathrm {i} \omega _\ell, \hat {\boldsymbol{q}}_{\ell } ]$ is determined as a solution of the following eigenvalue problem:

(2.6)\begin{equation} \boldsymbol{J}_{(\omega_\ell, m_\ell)} \hat{\boldsymbol{q}}_{(z_{\ell})} \equiv \left({\rm i}\omega_\ell \boldsymbol{B} - \left.\frac{\partial \boldsymbol{F}}{\partial \boldsymbol{q} }\right|_{\boldsymbol{q} = \boldsymbol{Q}_0, \boldsymbol{\eta} = \boldsymbol{0}} \right) \hat{\boldsymbol{q}}_{(z_{\ell})} = 0, \end{equation}

where the linear operator $\boldsymbol{J}$ is the Jacobian of (2.1), and $({\partial \boldsymbol{F} }/{\partial \boldsymbol{q} }|_{\boldsymbol{q} = \boldsymbol{Q}_0, \boldsymbol {\eta } = \boldsymbol {0}} ) \hat {\boldsymbol{q}}_{(z_{\ell })} = \boldsymbol{L}_{m_\ell } \hat {\boldsymbol{q}}_{(z_{\ell })} + \boldsymbol{N}_{m_\ell }(\boldsymbol{Q}_0, \hat {\boldsymbol{q}}_{(z_{\ell })} ) + \boldsymbol{N}_{m_\ell }( \hat {\boldsymbol{q}}_{(z_{\ell })} , \boldsymbol{Q}_0 )$. The subscript $m_\ell$ indicates the azimuthal wavenumber used for the evaluation of the operator. In the following, we account for eigenmodes $\hat {\boldsymbol{q}}_{(z_{\ell })} (r,z)$ that have been normalised in such a way $\langle \hat {\boldsymbol{u}}_{(z_\ell )}, \hat {\boldsymbol{u}}_{(z_\ell )} \rangle _{L^{2}} = 1$.

The identification of the core region of the self-excited instability mechanism (Gianneti & Luchini Reference Gianneti and Luchini2007) is evaluated by means of the structural sensitivity tensor

(2.7)\begin{equation} \boldsymbol{\textsf{S}}_s = ( \hat{\boldsymbol{u}}^{{\dagger}} )^\ast\otimes \hat{\boldsymbol{u}}. \end{equation}

2.2.2. Methodology for the study of mode selection

In the following, we briefly outline the main aspects of the methodology employed in the study of mode interaction or unfolding of a bifurcation with codimension-two, a comprehensive explanation is left to Appendix A. Herein, we use the concept of mode interaction as a synonym of the analysis of a bifurcation with codimension-two, that is, a bifurcation satisfying two conditions, e.g. a bifurcation where two modes become at the same time unstable. The determination of the attractor or coherent structure is explored within the framework of equivariant bifurcation theory. The trivial steady state is axisymmetric, i.e. the symmetry group is the orthogonal group $O(2)$. Near the onset of the instability, the dynamics can be reduced to that of the centre manifold. Particularly, due to the non-uniqueness of the manifold, one can always look for its simplest polynomial expression, which is known as the normal form of the bifurcation. The reduction to the normal form is carried out via a multiple scales expansion of the solution $\boldsymbol{Q}$ of (2.4). The expansion considers a two scale development of the original time $t \mapsto t + \varepsilon ^2 \tau$, here, $\varepsilon$ is the order of magnitude of the flow disturbances, assumed to be small $\varepsilon \ll 1$. In this study we carry out a normal form reduction via a weakly nonlinear expansion, where the small parameters are

(2.8a,b)\begin{equation} \varepsilon_{\delta_u}^2 = \delta_{u,c} - \delta_{u} \sim \varepsilon^2 \quad\text{and}\quad \varepsilon_{\nu}^2 = \left( \nu_c - \nu \right) = ( {{Re}}_c^{{-}1} - {{Re}}^{{-}1} ) \sim \varepsilon^2. \end{equation}

A fast time scale $t$ of the self-sustained instability and a slow time scale of the evolution of the amplitudes $z_i(\tau )$ are also considered in (2.13), for $i=1,2,3$. The ansatz of the expansion is as follows:

(2.9)\begin{equation} \boldsymbol{Q}(t,\tau) = \boldsymbol{Q}_0 + \varepsilon \boldsymbol{q}_{(\varepsilon)}(t,\tau) + \varepsilon^2 \boldsymbol{q}_{(\varepsilon^2)}(t,\tau) + O(\varepsilon^3). \end{equation}

Herein, we evaluate the mode interaction between two steady symmetry-breaking states with azimuthal wavenumbers $m_1 = 1$ and $m_2 = 2$, that is,

(2.10) \begin{align} \boldsymbol{q}_{(\varepsilon)}(t,\tau) & = ( z_1(\tau) \hat{\boldsymbol{q}}_{(z_{1})} (r,z) \exp({-\mathrm{i} m_1 \theta}) + \text{c.c.} ) \nonumber\\ &\quad + ( z_2(\tau) \hat{\boldsymbol{q}}_{(z_{2})} (r,z) \exp({-\mathrm{i} m_2 \theta}) + \text{c.c.} ), \end{align}

where $z_1$ and $z_2$ are the complex amplitudes of the two symmetric modes $\hat {\boldsymbol{q}}_{(z_{1})}$ and $\hat {\boldsymbol{q}}_{(z_{1})}$. Note that the expansion of the left-hand side of (2.4) up to third order is as follows:

(2.11)\begin{equation} \varepsilon \boldsymbol{B} \frac{\partial \boldsymbol{q}_{(\varepsilon)}} {\partial t} + \varepsilon^2 \boldsymbol{B} \frac{\partial \boldsymbol{q}_{(\varepsilon^2)}} {\partial t} + \varepsilon^3 \left[ \boldsymbol{B} \frac{\partial \boldsymbol{q}_{(\varepsilon^3)}} {\partial t}\right] + O(\varepsilon^4), \end{equation}

and the right-hand side, respectively,

(2.12)\begin{equation} \boldsymbol{F} (\boldsymbol{q}, \boldsymbol{\eta}) = \boldsymbol{F}_{(0)} + \varepsilon \boldsymbol{F}_{(\varepsilon)} + \varepsilon^2 \boldsymbol{F}_{(\varepsilon^2)} + \varepsilon^3 \boldsymbol{F}_{(\varepsilon^3)} + O(\varepsilon^4). \end{equation}

Then, the problem up to third order in $z_1$ and $z_2$ can be reduced to (Armbruster, Guckenheimer & Holmes Reference Armbruster, Guckenheimer and Holmes1988)

(2.13) \begin{equation} \left.\begin{gathered} \dot{z}_1 = \lambda_1 z_1 + e_3 \bar{z}_1 z_2 + z_1( c_{(1,1)} |z_1|^2 + c_{(1,2)} |z_2|^2 ), \\ \dot{z}_2 = \lambda_2 z_2 + e_4 z^2_1 + z_2( c_{(2,1)} |z_1|^2 + c_{(2,2)} |z_2|^2 ). \end{gathered}\right\} \end{equation}

where $\lambda _1$ and $\lambda _2$ are the unfolding parameters of the normal form. The procedure followed for the determination of the coefficients $c_{(i,j)}$ for $i,j=1,2$ and $e_3$ and $e_4$ is left to Appendix A. An exhaustive analysis of the nonlinear implications of this normal form on the dynamics is left to § 5.

4.1. Spectrum

Herein, we analyse the asymptotic linear stability of the steady state in four distinct configurations. The first spectrum, depicted in figure 8(a), has been computed for a velocity ratio $\delta _u = 1$. Similarly, the second spectrum corresponds to a velocity ratio $\delta _u=0.28$, which represents the middle branch after the saddle node, that is, the equivalent of the marker (b) in figure 4(b) for ${{Re}}=250$. These two configurations have been determined for a duct wall length $L=1$. The remaining two spectra have been computed for duct wall distances of $L=0.5$ and $L=2$, which are illustrated in figures 8(c) and 8(d), respectively. The computation of the spectrum reveals the existence of eigenmodes, with azimuthal wavenumbers $m=0$, $m=1$ and $m=2$, that become unstable.

Figure 8. Spectrum computed at four different configurations of (${{Re}},L,\delta _u$) for $m=0,1,2$. The inset inside each panel magnifies the region near the origin. Stationary or low frequency modes are designated $S$, while oscillating/flapping modes are designated $F$, with the azimuthal wavenumber as the subscript: (a) $Re=250$, $L=1$, $\delta _u=1$; (b) $Re=250$, $L=1$, $\delta _u=0.28$; (c) $Re=800$, $L=0.5$, $\delta _u=1$; (d) $Re=250$, $L=2$, $\delta _u=1$.

First, the four spectra display three types of continuous branches, referred to as $b_i$ ($i=1,2,3$), as was the case in the configuration of coaxial jets described by Canton et al. (Reference Canton, Auteri and Carini2017). The branch $b_3$ is composed of spurious modes. The branch $b_2$ is constituted of modes localised within the jet shear layers. While the branch $b_1$ is composed by nearly steady modes with support in the fluid region surrounding the jets.

Second, in the four configurations we find two non-oscillating unstable modes (or nearly neutral, as is the case in figure 8c) with azimuthal wavenumbers $m=1$ and $m=2$, hereinafter referred to as modes $S_1$ and $S_2$, respectively. These two modes are depicted in figure 9, which illustrates their axial velocity components for the four configurations. Their spatial distribution is mostly localised inside the recirculating region of the flow, but they are also supported along the shear layer of the jets. Evaluating both the direct and adjoint modes, we can identify the core of the global instability from the maximum values of the function $||{\boldsymbol{\mathsf{S}}}_s(r,z)||_F$, which has been defined in (2.7).

Figure 9. Axial velocity component of the non-oscillating global modes $S_1$ (bottom half of each panel) and $S_2$ (top half of each panel). The label of the panels coincides with the label of figure 8: (a) $Re=250$, $L=1$, $\delta _u=1$; (b) $Re=250$, $L=1$, $\delta _u=0.28$; (c) $Re=800$, $L=0.5$, $\delta _u=1$; (d) $Re=250$, $L=2$, $\delta _u=1$.

Figure 10 illustrates the sensitivity maps for the modes displayed in figure 9(a,c,d). The sensitivity maps $||{\boldsymbol{\mathsf{S}}}_s(r,z)||_F$ are compact supported within the region of recirculating fluid, featuring negligible values elsewhere. The maximum values of the sensitivity maps, displayed in figure 10(a,c,e) for the mode $S_1$, are found within the inner vortex ring, in particular near the downstream part of the inner vortical region, and on the interface between the two vortical rings. By increasing the wall length separating the jet streams, the wavemaker moves downstream towards the right end of the inner vortical region. A similar observation is drawn from figure 10(b,d, f), where the core of the instability is also found within the inner vortex ring. Similar observations were drawn in the case of the wake behind rotating spheres (Sierra-Ausín et al. Reference Sierra-Ausín, Lorite-Díez, Jiménez-González, Citro and Fabre2022), where the core of the instability was also found near the downstream part of the recirculating flow region. Therein, it was concluded that the instability is supported by the recirculating flow region.

Figure 10. Structural sensitivity map $||{\boldsymbol{\mathsf{S}}}_s(r,z)||_F$. White lines are employed to represent the steady-state streamlines: (a) $S_1$ ($Re=250, L=1, \delta _u=1$); (b) $S_2$ ($Re=250, L=1, \delta _u=1$); (c) $S_1$ ($Re=800, L=0.5, \delta _u=1$); (d) $S_2$ ($Re=800, L=0.5, \delta _u=1$); (e) $S_1$ ($Re=250, L=2, \delta _u=1$); ( f) $S_2$ ($Re=250, L=2, \delta _u=1$).

Figure 8(d) illustrates the existence of two oscillating/flapping modes with azimuthal wavenumbers $m=1$ and $m=2$, hereinafter referred to as $F_1$ and $F_2$, respectively. The axial velocity component of these two modes is displayed in figure 11, together with their associated structural sensitivity map. The unsteady modes $F_1$ and $F_2$ possess a much larger spatial support than $S_1$ and $S_2$. They are formed by an array of counter-rotating vortex spirals sustained along the shear layer of the base flow. For the mode $F_2$ the amplitude of these structures grows downstream of the nozzle, in the axial direction, with a maximum around $z\approx 70$, after which they slowly decay. The mode $F_1$ grows further downstream, with a maximum around $z \approx 300$. The spatial structure of these eigenmodes resembles the axisymmetric mode of figure 9 in Canton et al. (Reference Canton, Auteri and Carini2017) or the optimal response modes determined by Montagnani & Auteri (Reference Montagnani and Auteri2019). As was the case for the non-oscillating modes, the core of the instability is found near the downstream part of the inner vortex ring. Tentatively, one may conclude that vortical perturbations are produced within the recirculating flow region and convected downstream while experiencing a considerable spatial amplification, which in turn justifies the resemblance to the optimal response modes determined by Montagnani & Auteri (Reference Montagnani and Auteri2019).

Figure 11. Axial velocity component of the oscillating global modes $F_1$ (a) and $F_2$ (b). Structural sensitivity map $||{\boldsymbol{\mathsf{S}}}_s(r,z)||_F$ of mode $F_1$ (c) and $F_2$ (d). White lines are employed to represent the steady-state streamlines.

There is an unstable $m=0$ mode, hereinafter referred to as $S_0$, in the spectrum displayed in figure 8(b). Such a mode, which is illustrated figure 12(a), is the result of a saddle-node bifurcation leading to the existence of multiple steady states, a feature that has been discussed in § 3. It is a mode that promotes the formation of a recirculating flow region attached to the duct wall. In § 3 we have remarked that the $S_0$ modes can be related to changes in the topology of the flow, and to a downstream shift of the recirculation bubble. Thus, it is not surprising that the core of the instability, shown in figure 12(b), is found on the interface between the recirculating region attached to the wall and the large recirculation bubble, and mostly in a region close to the axis found near the leftmost end of the recirculation bubble. The changes in the base flow due to the $S_0$ mode have an impact on the instability core of the $S_1$ and $S_2$ modes, which are depicted in figures 12(c) and 12(d), respectively. The maximum values of the structural sensitivity are found on the leftmost end of the recirculation bubble near the axis of revolution, while it is found in the centre of the recirculation bubble for the mode $S_2$.

Figure 12. (a) Global mode $S_0$ for the configuration ($Re=250, L=1, \delta _u=0.28$). The top half of (a) represents the axial velocity, while the bottom half depicts the radial velocity component. Structural sensitivity map $||{\boldsymbol{\mathsf{S}}}_s(r,z)||_F$ of the modes $S_0$ (b), $S_1$ (c) and $S_2$ (d). White lines are employed to represent the steady-state streamlines.

4.2. Annular jet configuration $\delta _u=0$

Herein, we investigate the effect of the duct wall length ($0.5 < L < 4$) on the linear stability of the annular jet ($\delta _u=0$).

The linear stability findings are summarised in figure 13, which displays the evolution of the critical Reynolds number with respect to the duct wall distance ($L$) for the four most unstable modes: two non-oscillating $S_1$ and $S_2$, and two oscillating $F_1$ and $F_2$. A cross-section view at $z=1$ is displayed in figure 14. Please note that, for the chosen set of parameters, the axisymmetric unsteady mode $F_0$ is always found at larger Reynolds numbers than the aforementioned modes, that is why, in the following, we only include the results for the $S_1$, $S_2$, $F_1$ and $F_2$ modes. This is one of the major differences from the case studied by Canton et al. (Reference Canton, Auteri and Carini2017). For small values of the duct wall length ($L \approx 0.1$) separating the jet streams, the dominant instability is a vortex-shedding mode, which in our nomenclature is referred to as $F_0$. On the contrary, for duct wall lengths in the interval $0.5 < L < 4$, the primary instability of the annular jet is a steady symmetry-breaking bifurcation that leads to a jet flow with a single symmetry plane, displayed in figure 14(a). In contrast, bifurcations that lead to the mode $S_2$ possess two orthogonal symmetry planes, see figure 14(b). In § 4.1 it has been established that non-oscillating modes $S_1$ and $S_2$ for $\delta _u = 1$ display the highest intensity within the region of recirculating fluid. Likewise, in the annular jet configuration, figure 15 demonstrates that the spatial distribution of these two stationary modes $S_1$ and $S_2$ is found inside the recirculation bubble. For jet distances $L < 2$, the second mode that bifurcates is $F_1$ mode, depicted in figure 16(a). This situation corresponds to a bifurcation scenario similar to other axisymmetric flow configurations, such as the flow past a sphere or a disk (Meliga, Chomaz & Sipp Reference Meliga, Chomaz and Sipp2009; Auguste, Fabre & Magnaudet Reference Auguste, Fabre and Magnaudet2010). For larger distances between jets, the scenario changes. The second bifurcation from the axisymmetric steady state is the $F_2$, displayed in figure 16(b). Other configurations where the primary or secondary instability involves modes with azimuthal component $m=2$ are swirling jets (Meliga, Gallaire & Chomaz Reference Meliga, Gallaire and Chomaz2012) and the wake flow past a rotating sphere (Sierra-Ausín et al. Reference Sierra-Ausín, Lorite-Díez, Jiménez-González, Citro and Fabre2022). The unsteady modes $F_1$ and $F_2$ display a similar structure to the unsteady modes discussed in § 4.1. They are formed by an array of counter-rotating vortex spirals developing in the wake of the separating duct wall and convected downstream, while experiencing an important spatial amplification until they eventually decay after reaching a maximum amplitude.

Figure 13. Linear stability boundaries for the annular jet ($\delta _u = 0$). (b) Frequency evolution of the unsteady modes. Legend: $S_1$ mode is displayed with a solid black line, $S_2$ with a solid red line and the $F_1$ and $F_2$ modes are depicted with dashed black and red lines, respectively.

Figure 14. Cross-sectional view at $z=1$ of the four unstable modes at criticality for the annular jet case ($\delta _u=0$). The streamwise component of the vorticity vector $\varpi _z$ is visualised by colour. (a) Mode $S_1$ for $L=0.5$. (b) Mode $S_2$ for $L=0.5$. (c) Mode $F_1$ for $L=3$. (d) Mode $F_2$ for $L=3$.

Figure 15. Global modes $S_1$ (a) and $S_2$ (b) at criticality for $L=0.5$ and $\delta _u = 0$. The top halves of each panel represent the axial velocity, while the bottom halves depict the radial velocity component. Black lines represent the streamlines of the base flow.

Figure 16. Axial velocity component of the neutral modes for $L=3$ and $\delta _u=0$; (a) $F_1$, (b) $F_2$.

4.3. Fixed distance between jets and variable velocity ratio $\delta _u$

In the following, we focus on the influence of the velocity ratio $\delta _u$ between jets for fixed jet distances $L$. Figure 17 displays the neutral curve of stability for jet distances (a) $L=0.5$ and (b) $L=1$. One may observe that the primary bifurcation is not always associated with the mode $S_1$ as is the case for $\delta _u=0$. For sufficiently large velocity ratios, the primary instability leads to a non-axisymmetric steady state with a double helix, corresponding to the unstable mode $S_2$. As can be appreciated in figure 9(b), for small values of $\delta _u$, the mode $S_1$ expands downstream over a relatively large area, having a higher activity than mode $S_2$, which is confined to the recirculation region. As the ratio between velocities is increased, as observed in figure 9(a), mode $S_2$ enlarges and resembles mode $S_1$, controlling the instability mechanism for large values of $\delta _u$. Another interesting feature, which could motivate a control strategy, is the occurrence of vertical asymptotes. This sudden change in the critical Reynolds number is due to the retraction, disappearance of the recirculation bubble and the formation of a new recirculating flow region, aspects that have been covered in § 3. For $L=0.5$, this sudden change occurs for $\delta _u \approx 0.25$, and for higher values of $\delta _u$ the critical Reynolds number is around twice larger than the one of the annular jet ($\delta _u=0$). The case of jet distance $L=1$ was discussed in § 3. The sudden change in the stability of the branch $S_1$ occurs in the range $\delta _u \in [0.25,0.5]$. Within this narrow interval, the primary branch of instability is the $F_1$. At around $\delta _u=0.4$, the primary bifurcation is again the branch $S_1$, which becomes secondary at around $\delta _u \approx 0.8$ in favour of the branch $S_2$. In figure 17 we have highlighted the codimension two point interaction between the $S_1-S_2$ modes, which will be analysed in detail in § 5. Around this point, we can observe the largest ratio (${Re_c|_{\delta _u\neq 0}}/{Re_c|_{\delta _u=0}}$) between the value of the critical Reynolds number of the primary instability for a concentric jet configuration ($\delta _u \neq 0$) and the annular jet problem ($\delta _u=0$).

Figure 17. Linear stability boundaries for the concentric jets (a) $L = 0.5$ and (b) $L = 1$. Same legend as figure 13.

4.4. Fixed velocity ratio $\delta _u$ and variable distance between jets

Figure 18 compares the results obtained for a constant velocity ratio when varying the distance between jets. As observed before, the increase of the distance between the jets has a de-stabilising effect. The largest critical Reynolds number is found at $\delta _u=0$, and the critical Reynolds number decreases with the duct wall length $L$ between the jet streams. The points of codimension two are highlighted in figure 18. We can appreciate that the interaction between the branch $S_1$ and $S_2$ happens for every velocity ratio $\delta _u$ explored, and it is the mode interaction associated with the smallest distance between jets. Additionally, for a velocity ratio $\delta _u=0.5$ there exist two points where the branches of the linear modes $S_1$ and $F_1$ intersect. Another feature of the neutral curves is the existence of turning points, which are associated with the existence of saddle-node bifurcations of the axisymmetric steady state, addressed in § 3. The saddle-node bifurcations of the steady state induce the existence of regions in the neutral curves with a tongue shape. These saddle-node bifurcations are also responsible for the formation of the vertical asymptotes observed in figure 17. Finally, it is of interest the transition of the modes $S_1$ and $S_2$, which induce the symmetry breaking of the axisymmetric steady state to slow low frequency spiralling structures. These modes have been identified for $\delta _u=0.5$ for $m=1$, $\delta _u=1$ for $m=2$ and $\delta _u=2$ for both $m=1$ and $m=2$. As it will be clarified in § 5, these oscillations issue from the nonlinear interaction of modes, emerging simultaneously for a specific Reynolds number, and changing their position as the most unstable global mode of the flow.

Figure 18. Neutral lines of the four modes found by studying the configuration of two concentric jets and fixing the velocity ratio; (a,b) $\delta _u =0.5$, (c,d) $\delta _u = 1$, (e, f) $\delta _u = 2$. Black lines: modes with $m=1$, red lines: modes with $m=2$. Straight lines: steady modes, dashed lines: unsteady modes. The discontinuity points, i.e. the points where the second most unstable mode (of a given type) becomes the most unstable, are highlighted with square markers.

5.1. Normal form, basic solutions and their properties

The linear diagrams of § 4 have shown the existence of the mode interaction between the modes $S_1$ and $S_2$. It corresponds roughly to the mode interaction that occurs at the largest critical Reynolds number for any value of $L$ herein explored. In this section, we analyse the dynamics near the $S_1:S_2$ organising centre. We perform a normal form reduction, which allows us to predict non-axisymmetric steady, periodic, quasiperiodic and heteroclinic cycles between non-axisymmetric states.

The mode interaction that is herein analysed corresponds to a steady–steady bifurcation with $O(2)$ symmetry and with strong resonance 1 : 2. Such a bifurcation scenario has been extensively studied in the past by Dangelmayr (Reference Dangelmayr1986), Jones & Proctor (Reference Jones and Proctor1987), Porter & Knobloch (Reference Porter and Knobloch2001), Armbruster et al. (Reference Armbruster, Guckenheimer and Holmes1988) and for the reflection symmetry-breaking case (SO(2)) by Porter & Knobloch (Reference Porter and Knobloch2005). In order to unravel the existence and the stability of the nonlinear states near the codimension two point, let us write the flow field as

(5.1) \begin{equation} \boldsymbol{q} = \boldsymbol{Q}_0 + {Re} [ r_1(\tau) {\rm e}^{{\rm i} \phi_1(\tau)} {\rm e}^{-{\rm i}\theta} \hat{\boldsymbol{q}}_{s,1} ] + {Re} [ r_2(\tau) {\rm e}^{{\rm i} \phi_2(\tau)} {\rm e}^{{-}2{\rm i} \theta} \hat{\boldsymbol{q}}_{s,2} ], \end{equation}

in polar coordinates for the complex amplitudes $z_1 = r_1 \textrm {e}^{\mathrm {i} \phi _1}$ and $z_2 = r_2 \textrm {e}^{\mathrm {i} \phi _2}$ where $r_j$ and $\phi _j$ for $j=1,2$ are the amplitude and phase of the symmetry-breaking modes $m=1$ and $m=2$, respectively. The complex-amplitude normal form (2.13) is expressed in this reduced polar notation as follows:

(5.2a)$$\begin{gather} \dot{r}_{1} = e_3 r_{1} r_{2} \cos(\chi) + r_{1} ( \lambda_{(s,1)} + c_{(1,1)} r_{1}^2 + c_{(1,2)} r_{2}^2 ), \end{gather}$$

(5.2b)$$\begin{gather}\dot{r}_{2} = e_4 r_{1}^2 \cos(\chi) + r_{2} ( \lambda_{(s,2)} + c_{(2,1)} r_{1}^2 + c_{(2,2)} r_{2}^2 ), \end{gather}$$

(5.2c)$$\begin{gather}\dot{\chi} ={-} \left( 2 e_3 r_{2} + e_4 \frac{r_{1}^2}{r_{2}} \right) \sin(\chi), \end{gather}$$

where the phase $\chi = \phi _2 - 2\phi _1$ is coupled to the amplitudes $r_1$ and $r_2$ because of the existence of the 1 : 2 resonance. The individual phases evolve as

(5.3)\begin{equation} \left.\begin{gathered} \dot{\phi}_1 = e_3 r_2 \sin(\chi), \\ \dot{\phi}_2 ={-} e_4 \frac{r_1^2}{r_2} \sin(\chi). \end{gathered}\right\} \end{equation}

Before proceeding to the analysis of the basic solutions of (5.2), we can simplify these equations by the rescaling

(5.4)\begin{equation} \left( \frac{r_1}{|e_3 e_4|^{1/2}}, \frac{r_2}{e_3} \right) \to (r_1,r_2), \end{equation}

which yields the following equivalent system:

(5.5a)$$\begin{gather} \dot{r}_{1} = r_{1} r_{2} \cos(\chi) + r_{1} ( \lambda_{(s,1)} + c_{11} r_{1}^2 + c_{12} r_{2}^2 ), \end{gather}$$

(5.5b)$$\begin{gather}\dot{r}_{2} = s r_{1}^2 \cos(\chi) + r_{2} ( \lambda_{(s,2)} + c_{21} r_{1}^2 + c_{22} r_{2}^2 ), \end{gather}$$

(5.5c)$$\begin{gather}\dot{\chi} ={-} \frac{1}{r_2} ( 2 r_{2}^2 + s r_1^2) \sin(\chi), \end{gather}$$

where the coefficients

(5.6a–e)\begin{align} s = \text{sign}(e_3 e_4), \quad c_{11} = \frac{c_{(1,1)}}{|e_3 e_4|}, \quad c_{12} = \frac{c_{(1,2)}}{e_3^2}, \quad c_{21} = \frac{c_{(2,1)}}{|e_3 e_4|}, \quad c_{22} = \frac{c_{(2,2)}}{e_3^2}. \end{align}

Finally, we consider a third normal form equivalent to the previous ones but which removes the singularity of (5.2) and (5.5) when $r_2 = 0$. Standing waves ($\sin \chi = 0$) naturally encounter this type of artificial singularity, which manifests in (5.5) as an instantaneous jump from one standing subspace to the other by a ${\rm \pi}$-translation. This is the case of the heteroclinic cycles, previously studied by Armbruster et al. (Reference Armbruster, Guckenheimer and Holmes1988) and Porter & Knobloch (Reference Porter and Knobloch2001). The third normal form, which we shall refer to as the reduced Cartesian normal form, takes advantage of the simple transformation $x=r_2 \cos (\chi )$, $y=r_2 \sin (\chi )$ (Porter & Knobloch Reference Porter and Knobloch2005)

(5.7a)$$\begin{gather} \dot{r}_{1} = r_{1} ( \lambda_{(s,1)} + c_{11} r_{1}^2 + c_{12} (x^2+y^2) + x ), \end{gather}$$

(5.7b)$$\begin{gather}\dot{x} = s r_{1}^2 + 2 y^2 + x( \lambda_{(s,2)} + c_{21} r_{1}^2 + c_{22} (x^2+y^2) ), \end{gather}$$

(5.7c)$$\begin{gather}\dot{y} ={-} 2 xy + y ( \lambda_{(s,2)} + c_{21} r_{1}^2 + c_{22} (x^2+y^2) ). \end{gather}$$

In this final representation, standing wave solutions are contained within the invariant plane $y=0$, and due to the invariance of (5.7) under the reflection $y \mapsto -y$, one can restrict attention, without loss of generality, to solutions with $y \geq 0$, cf. Porter & Knobloch (Reference Porter and Knobloch2001).

The system (5.5) possess four types of fixed points, which are listed in table 1.

Table 1. Definition of the fixed points of the reduced polar normal form (5.5); $\sigma _{\pm }$ is defined in (5.8), the polynomial $P_{MM}$ is defined in (5.10) and $\varSigma _{TW} \equiv 4 c_{11} + 2 (c_{12} + c_{21}) + c_{22}$.

First, the axisymmetric steady state (O) is represented by $(r_1,r_2) = (0,0)$, so it is the trivial steady state of the normal form. The second steady state is what it is denoted as the pure mode (P). In the original coordinates, it corresponds to the symmetry-breaking structure associated with the mode $S_2$. This state bifurcates from the axisymmetric steady state (O) when $\lambda _{(s,2)} = 0$. The third fixed point is the mixed-mode state (MM), which is listed in table 1. It corresponds to the reflection symmetry preserving state associated with the mode $S_1$. It may bifurcate directly from the trivial steady state O, when $\lambda _{(s,1)} = 0$ or from P whenever $\sigma _{+} = 0$ or $\sigma _{-} = 0$, where $\sigma _{\pm }$ is defined as

(5.8)\begin{equation} \sigma_{{\pm}} \equiv \lambda_{(s,1)} - \frac{-\lambda_{(s,2)} c_{12}}{c_{22}} \pm \sqrt{\frac{-\lambda_{(s,2)}}{c_{22}}} . \end{equation}

The representation in the reduced polar form is

(5.9a,b)\begin{equation} r_{1,MM} ={-} \frac{ \lambda_{(s,1)} \pm r_{2,MM} + c_{12} r_{2,MM}^2 }{c_{11}}, \quad \cos(\chi_{MM}) ={\pm} 1, \end{equation}

and the condition ${P}_{{MM}}(r_{2,MM} \cos (\chi _{MM})) = 0$, where ${P}_{{MM}}$ is defined as

(5.10)\begin{align} {P}_{{MM}}(x) \equiv s \mu_1 + (s + c_{21} \lambda_{(s,1)} - c_{(1,1)} \lambda_{(s,2)}) x + (c_{21} + s c_{12}) x^2 + (c_{12}c_{21} - c_{11}c_{22})x^3. \end{align}

Finally, the fourth fixed point of the system represents travelling waves (TW). It is surprising that the interaction between two steady states causes a time-periodic solution. The travelling wave emerges from MM in a parity-breaking pitchfork bifurcation that breaks the reflection symmetry when $\cos (\chi _{TW}) = \pm 1$. The TW drifts at a steady rotation rate $\omega _{TW}$ along the group orbit, i.e. the phases $\dot {\phi }_1 = r_{2,TW} \sin (\chi _{TW})$ and $\dot {\phi }_2 = -s ({r_{1,TW}^2}/{r_{2,TW}}) \sin (\chi _{TW})$ are non-null.

Mixed modes and travelling waves may further bifurcate into standing waves (SW) and modulated travelling waves (MTW), respectively. These are generic features of the 1 : 2 resonance for small values of $\lambda _{(s,1)}$ and $\lambda _{(s,2)}$, when $s = -1$. In the original coordinates, SW are periodic solutions, whereas MTW are quasiperiodic. Standing waves emerge via a Hopf bifurcation from MM when the conditions ${P}_{{SW}} (r_{2,MM} \cos (\chi _{MM}) ) > 0$ for

(5.11)\begin{equation} {P}_{{SW}} (x) \equiv (2 c_{22} x^3 - s r_1^2) c_{11} - (2 c_{12} x + 1)(c_{21}x + s)x, \end{equation}

and the one listed in table 2 are satisfied. The MTW are created when a torus bifurcation happens on the travelling wave branch when the conditions listed in table 2 are satisfied.

Table 2. Definition of the limit cycles of the reduced polar normal form (5.5).

Another remarkable feature of (5.2) is the existence of robust heteroclinic cycles that are asymptotically stable. When $s = -1$, there are open sets of parameters where the reduced polar normal form exhibits structurally stable connections between ${\rm \pi}$-translations on the circle of pure modes, cf. Armbruster et al. (Reference Armbruster, Guckenheimer and Holmes1988). These structures are robust and have been observed in a large variety of systems, (Palacios et al. Reference Palacios, Gunaratne, Gorman and Robbins1997; Mercader, Prat & Knobloch Reference Mercader, Prat and Knobloch2002; Nore et al. Reference Nore, Tuckerman, Daube and Xin2003; Mariano & Stazi Reference Mariano and Stazi2005; Nore, Moisy & Quartier Reference Nore, Moisy and Quartier2005). In addition to these robust heteroclinic cycles connecting pure modes, there exist more complex limit cycles connecting O, P, MM and SW, cf. Porter & Knobloch (Reference Porter and Knobloch2001). These cycles are located for larger values of $\lambda _{(s,1)}$ and $\lambda _{(s,2)}$, with a possibly chaotic dynamics (Shilnikov type). In this study, we have not identified any of these. Finally, a summary of the basic solutions and the bifurcation path is sketched in figure 19.

Figure 19. Schematic representation of the basic solutions of (5.2) and their bifurcation path.