A Numerical Study of the Global Instability in Counter-Current Homogeneous Density Incompressible Round Jets

The paper focuses on a global instability phenomenon in counter-current round jets issuing from co-axial nozzles. Three different configurations that differ in a way of the counter flow generation are investigated. Besides typical configurations used in experimental and numerical research performed so far, in which suction applied in an annular nozzle is a driving force for the counterflow, a novel set-up is proposed where the annular nozzle is oriented in the opposite direction and placed above the main one. Such a configuration eliminates the suction of fluid from the main jet, which in previous research was found to have a destructive impact on the occurrence of the global flow instability. The research is performed using a large eddy simulation (LES) method and the computations are carried out applying a high-order numerical code, the accuracy of which has been proven in previous works and also in the present research through comparisons with available experimental data. The research is complemented by the linear stability analysis which supports the LES results and formulated conclusions. In agreement with a number of the previous works it has been shown that the global modes can be triggered only when the velocity ratio (I) between the main jet velocity and the velocity of the jet issuing from the annular nozzle is above a certain threshold level (Icr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{\text {cr}}$$\end{document}). It has been shown that in the classical configurations of the co-axial nozzles the range of I≥Icr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I\ge I_{\text {cr}}$$\end{document} for which the global instability phenomenon exists is very narrow and it disappears for larger velocity ratios. Reasons for that have been identified through detailed scrutiny of instantaneous flow pictures. In the new set-up of the nozzles the global instability persists for a significantly wider range of I. It has been shown that Icr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{\text {cr}}$$\end{document} depends on both the momentum thickness of the mixing layer formed between the counter-current streams and the applied configuration of the nozzles. The LES results univocally showed that the latter factor decides on the type of the instability mode (Mode I or Mode II) that emerges in the flow, as it directly influences on a length of the region where the counter-current streams are parallel allowing the growth of short or long wave disturbances characteristic for Mode I and Mode II, respectively.


Introduction
Round jets, due to their occurrence in many industrial applications, have received huge interest for decades. Special attention has been devoted to an enhancement of the mixing processes, which are essential for an efficient and safe work of all devices where proper mixing is crucial for a good performance, e.g. in combustion processes where the jets issue from fuel nozzles, heat exchangers where the jets are used for surface cooling/heating, boundary layer control devices with the jets acting as vortex generators, or in various types of mixers and injectors. The dynamics of the flow in the round jets seem to be well recognised and discussed in details in review papers (Ball et al. 2012;Lipari and Stansby 2011). The understanding of the transition mechanisms and flow behaviour in the round jets facilitates sophisticated research on their control involving artificial intelligence techniques , an extremum seeking method  or predetermined approaches (Gohil et al. 2015;Tyliszczak 2018).
An interesting phenomenon emerging in the round jets, which significantly intensifies the mixing intensity, is the self-excited global instability phenomenon triggered by absolutely unstable local flow regions (Huerre and Monkewitz 1990). Throughout the paper, an absolute mode or absolute instability will be understood as the instability of a parallel base flow for which the spatio-temporal stability theory shows the positive temporal growth rate. Flow oscillations, triggered by a locally absolutely unstable flow region leading to the selfsustained mechanism, will be understood as global modes and global oscillations. Since the 1980s the global instability phenomenon has been studied in several investigations in various configurations using both experimental and numerical methods. Although literature in this field is very extensive, see for instance (Monkewitz and Sohn 1988;Strykowski and Niccum 1991;Liang and Maxworthy 2005;Lesshafft et al. 2006;Oberleithner et al. 2011;Meliga et al. 2012;Tammisola and Juniper 2016;Wawrzak et al. 2019;Vanierschot et al. 2020), the global instability phenomenon and the flow behaviour it induces remains not fully recognised.
The linear stability calculations (Monkewitz and Sohn 1988;Jendoubi and Strykowski 1994) showed that the absolute instability can be observed in variable density jets and constant density counter-current jet configurations. Concerning the variable density jets extensive experimental, numerical and theoretical investigations have been done taking into account the impact of operating conditions (Reynolds number, momentum thickness and jet-ambient fluid density ratio) (Monkewitz and Sohn 1988;Monkewitz et al. 1990; Kyle and Sreenivasan 1993;Hallberg and Strykowski 2006;Lesshafft et al. 2006;Lesshafft and Marquet 2010;Srinivasan and Strykowski 2010;Coenen and Sevilla 2012;Foysi et al. 2010;Boguslawski et al. 2016;Zhu et al. 2017), for the last decades. Global modes have been shown to exist over a wide range of these parameters, but the issue of universality of the flow response has not been fully discovered. The problem of global oscillations in constant density counter-current configurations, however, remains the least explored. Jendoubi and Strykowski (1994) demonstrated, applying the linear stability theory (LST), that an external counterflow imposed on the jet triggered the absolute instability phenomenon if the velocity of the counter-flowing jet was higher than 14% of the main jet. The two distinct absolutely unstable modes called Mode I and Mode II that occur in a certain range of the control parameters were identified. Both modes were found to become more unstable as the velocity ratio was increased, the density ratio was decreased or the Mach number was decreased. Besides, Jendoubi and Strykowski (1994) showed that the pressure fluctuations of Mode I has a maximum in the shear layer and vanishes at the centreline whereas in the case of Mode II the perturbation achieves a maximum at the jet centreline and decays monotonically in the radial direction. The stability problem has subsequently been explored by a number of researchers. Increasing CPU capability, as well as the development of appropriate numerical techniques, allowed scientists to avoid a priori assumptions of a quasi-parallel base state and conduct a global linear stability analysis. Concerning the jet configurations, several global linear stability analyses have been performed to study the critical parameters beyond which the onset of the global instability appears (Nichols and Lele 2011;Garnaud et al. 2013a, b;Coenen et al. 2017;Qadri et al. 2018) yielding predictions in fair agreement with experimental observations. Among others, the linear frequency response analysis performed by Coenen et al. (2017) revealed that the small perturbations may suffer strong amplifications producing the flow pattern similar to a globally unstable jet. Hypothetically, in an experiment the low level noise could therefore be amplified and observed as sustained coherent wavepackets (Coenen et al. 2017). To the best authors knowledge, the newer stability techniques have not been applied to study counter-current jet configurations. Despite their application lie out of the present work they use in the future may shed light on the results discussed in the paper.
Up to date, in all the experimental research the counter-current round jets were generated using two co-axial nozzles by placing an annular nozzle around the main jet nozzle and sucking the fluid from the region surrounding a jet. Strykowski and Niccum (1991) demonstrated a great potential of such a counterflow configuration in terms of the mixing enhancement. They analysed the counter-current round jets for a wide range of the velocity ratio I = −U 2 ∕U 1 = 0.0 − 0.2 ( U 1 -the jet velocity, U 2 -the velocity of the counterflow). If the velocity ratio was above a critical value I cr ≈ 0.14 well-organised vortices were formed at a single characteristic frequency, which was the most pronounced in the shear layer. As a consequence of the organised vortex motion, the level of the velocity fluctuations, as well as the decay rate of the mean velocity at the jet axis, were reduced. They argued that the transition at I cr was a result of the self-excited global oscillations that originated from the absolute-convective instability transition. Jendoubi and Strykowski (1994) compared the results of the stability analysis and experimental data of Strykowski and Niccum (1991) and suggested that the observed oscillations originated from the absolute Mode I. Strykowski and Wilcoxon (1993) performed an analogous experiment, however, they used an extension tube put on the external nozzle to increase a distance along which the main and counter-current jet could interact through a shear layer formed between them. It turned out that this modification resulted in a more pronounced effect on the overall jet development and global mixing enhancement of the flow. They also presented the profiles of time-averaged velocity fluctuations, which were significantly different from those shown by Strykowski and Niccum (1991). The fluctuations levels were elevated and their maxima were shifted upstream when the velocity ratio was above I cr .
Some light on the dependence of the jet behaviour on the initial conditions and discrepancies between the results reported by Strykowski and Niccum (1991) and Strykowski and Wilcoxon (1993) was shed by the linear stability analysis. The research of Jendoubi and Strykowski (1994) showed that in the experiment of Strykowski and Niccum (1991) the absolute Mode I was identified as the global oscillations were observed mainly in the shear layer. Contrary, Strykowski and Wilcoxon (1993) observed the disturbances throughout the entire jet region that were related to the appearance of the global Mode II. Later experimental research of Favre-Marinet and , Boguslawski et al. (1999) and Asendrych (2007) revealed the formation of side-jets and the substantial increase of the jet spreading rate, which were commonly regarded as a manifestation of the global instability regimes.
Although the research cited above univocally showed that the counterflow induces the global instability phenomenon, none of these works unambiguously identified its origin and a few questions without precise answers are left. The most important are: Why the global instability that occurs for I ≥ I cr disappears for large velocity ratios, as observed by Boguslawski et al. (1999)? Why adding the extension tube changes the mode of the global oscillations from Mode I to Mode II? Can we effectively control the appearance of these modes and extend the range of the global oscillations for large velocity ratios? The present work attempts to dispel these doubts.
We apply the LES method, which similarly to the direct numerical simulation (DNS), seems to be an ideal tool for such a study. LES and DNS turned out to be very useful for predictions of the global instability in low-density jets Foysi et al. 2010;Boguslawski et al. 2016), annular jets ), counter-current mixing layers formed in square and rectangular jets (Grinstein and DeVore 2002) and planar mixing layers (Kalghatgi et al. 2011;Kalghatgi and Acharya 2012), as well as in advanced studies on dependence of jet flows on initial conditions (Tyliszczak and Geurts 2019) and their control, e.g. in research on bifurcating (Tyliszczak and Geurts 2014) and blooming jets (Tyliszczak 2015;Gohil et al. 2015;Tyliszczak 2018). A variety of these researches and thorough comparisons with experimental data show that LES and DNS can be treated as reliable approaches also for studies focused on the global instability phenomenon in counter flowing round jets complexity of which does not differ significantly from the cases analysed in the studies cited above. In the paper, the evidence for that will be comparisons with available measurements and results from the linear stability theory.
The paper is organised as follows. The next section describes the nozzle geometries, provides details of the computational configuration, boundary conditions, and applied numerical method. Section 3 shows the results of the linear stability calculations to which we will refer throughout the paper, and the LES results are presented in Sect. 4. Conclusions and outlook for further research directions are formulated in Sect. 5.

Computational Configurations and Numerical Method
In this paper, we study the jet flows in three different configurations of the nozzles labelled as CONF-1, CONF-2 and CONF-3 schematic drawings of which are shown in Fig. 1 along with the view of the computational domain. The jets issuing from the nozzles evolve in a cylindrical domain with dimensions L r and L x in the radial (r) and axial (x) directions, respectively. In CONF-1 and CONF-2 the annular nozzles are placed around the central main jet nozzle. These set-ups correspond to the experimental research of Strykowski and Niccum (1991) and Strykowski and Wilcoxon (1993), respectively. Suction of the flow through the co-axial nozzles results in the occurrence of the counter-flow region in the direct vicinity of the main jet. It causes that the fluid flowing to the suction slot originates not only from the jet surroundings but also from its periphery, as can be deduced from Fig. 2 showing instantaneous flow pictures in the direct vicinity of the nozzles. It turned out that for I > 0.06 the flow rate originated from the main jet increased rapidly , which caused a destruction of the region of the parallel flow that in LST is assumed as the base flow. Thus, the experimental results for larger I values cannot be compared with the theoretical predictions. To overcome this problem in the present studies we include the third configuration (CONF-3) in which the suction of the fluid from the main jet stream is eliminated. In CONF-3 the counter-current flow is generated by placing the annular nozzle above the main one and orient it in the opposite direction. Similar configurations are used in innovative MILD-type combustion chambers (Kruse et al. 2015). In such devices, the counter-flowing jet of flue gases creates a shear layer between the fresh fuel and oxidiser streams. The shear layer triggers the mixing of the exhaust gases and incoming reactants, which uniform dilution is a key parameter in the MILD-type combustion process. As can be seen in Fig. 3, in CONF-3 the fluid leaving the upper nozzle engulfs the main jet, forms the counter-current region, and smoothly flows radially on the inclined wall of the central nozzle. In a real MILD-type combustion chamber with the fuel issuing from the main jet and the oxidiser flowing outside, the global oscillations generated by the counter-flow could significantly enhance their mixing. Hence, the configuration CONF-3 is important both from the point of view of the fundamental studies but also for practical applications. The main flow parameters and characteristic dimensions of the nozzles are shown in Figs. 2 and 3. The computational domain is the cylinder with sizes L r = 6D 1∕2 and L x = 15D 1∕2 where D 1∕2 is the main nozzle diameter equal to D 1∕2 = 2R 1∕2 , where R 1∕2 is a radius at which the axial velocity is equal to the half velocity defined as U 1∕2 = (U 1 + U 2 )∕2 . Preliminary simulations showed that the size of the computational domain is sufficiently large, i.e., the impact of its boundaries on the jet dynamics is negligible. As shown in Fig. 1 the main nozzles are placed below the lower side of the cylinder and in the simulations, these nozzles are not included in the computational domain. Instead, in CONF-1 and in CONF-2 at the lower side of the cylinder the inlet velocity is prescribed by the hyperbolic tangent (HT) profile, which accurately reflects the real velocity distribution formed in a jet in a short distance after leaving a nozzle. In CONF-2 an extension tube, which the length and outer diameter are equal to D 1∕2 and 1.5D 1∕2 , is embedded in the computational domain and its presence is modelled using the immersed boundary method (Mittal and Iaccarino 2005). It allows mimicking the impact of complex shape solid objects on the flow in simulations performed on structured meshes. The main and opposite nozzles in CONF-3 are treated in the same way. In this case, rather than assuming the inlet velocity through the HT profile, we define the inlet velocity using the Blasius profile, which constitutes real nozzles with a high contraction. In CONF-3 the main parameters are the Fig. 3 Contours of the instantaneous axial velocity distribution in the main cross-section in the configuration CONF-3 dimension of the lower nozzle D and the velocities at the inlets to the main and opposite nozzles ( U IP , U IA ) relations of which with the flow parameters in CONF-1 and CONF-2 will be discussed in the following subsections.

LES Solver
We consider an incompressible homogeneous density flow for which the continuity and Navier-Stokes equations in the framework of LES are given as: where the overbar represents spatial filtering. The variables in the equations are the velocity vector , the density , the pressure P and the viscous stress tensor . The sub-filter stress tensor is given by sgs = 2 sgs , where is the rate of strain tensor of the resolved velocity field and sgs is the sub-grid viscosity which in the present work is modelled as proposed by Vreman (2004). The LES code employed in this study is called SAILOR. It is an academic high-order solver, which was previously verified in computations of isothermal, non-isothermal, excited and reacting jets (Boguslawski et al. 2017). The solution algorithm is formulated in the cylindrical ( x, r, ) coordinate system using a half-staggered pattern of pressure and velocity nodes distribution (Tyliszczak 2014(Tyliszczak , 2016. The solution algorithm is based on the projection method for pressure-velocity coupling. The time integration is based on a low storage 3rd order Runge-Kutta method (Williamson 1980). The spatial discretisation is performed applying a sixth order compact difference scheme (Lele 1992) for the radial and axial direction and the Fourier pseudospectral method (Canuto et al. 1988) for the azimuthal one. The immersed boundary method used to model the solid elements of the nozzle decreases the solution accuracy (Mittal and Iaccarino 2005) and its combination with high-order schemes formally results only in 2nd order accuracy (Tyliszczak and Ksiezyk 2018) due to local errors in the vicinities of the walls. However, it is expected that further from the walls the resolving efficiency of the high-order schemes is recovered and definitely better than that obtained using low-order schemes.

Boundary Conditions
Instantaneous inlet velocity profiles (u(r, t)) are defined through the supperposition of the mean velocity profiles ( U m (r) ) and fluctuating velocity components ( u turb (r, t) ) as u(r, t) = U m (r) + u turb (r, t) . The velocity fluctuations are computed applying the digital filtering method proposed by Klein et al. (2003), which guarantees a properly correlated turbulent velocity field. In CONF-1 the inlet velocity profile is defined by the hyperbolic tangent as: where R = (U 1 − U 2 )∕(U 1 + U 2 ) , is the momentum thickness and U c is the co-flowing velocity. The co-flow is added to mimic the natural entrainment from the surroundings and we assume U c = 0.05U 1 . It was shown in Da Silva and Métais (2002) that the co-flow at this level does not alter the jet dynamics. In the case of CONF-2 the inlet boundary includes the lower side of the extension tube and in this case, the mean velocity profile is given as: Figure 4 shows the mean inlet velocity profiles used in CONF-1 and CONF-2. The symbol denotes the distribution of the points in the computational mesh details of which are discussed later. The Blasius profile is defined based on the centreline velocity U IP and is obtained by solving the Blasius equation by the power series method which leads to the following formula: where The Blasius profile is applied in CONF-3 and also in selected cases in CONF-2 in order to match the inlet velocity distribution reported in the experimental research of Favre-Marinet and  to which the LES results are compared. In the upper part of the annular nozzle (CONF-3) a uniform velocity profile is assumed with the velocity U AI . Due to contraction of the nozzle, the flow accelerates and at the exit, its bulk velocity equals U AE = 7∕3U AI . The mean velocity profiles used in CONF-3 are shown in Fig. 5. In all the configurations at the side boundary, the streamwise velocity is taken equal to the co-flow velocity while the radial and azimuthal velocities are set equal to zero. The pressure at the inlet and side boundaries is computed using the Neumann condition ⋅ ∇p = 0 with being the outward normal vector. At the outlet plane, all the velocity components are computed from the convective boundary condition u i ∕ t + C u i ∕ x = 0 with C being the convection velocity which is calculated as the mean velocity in the outlet plane, at every time step. The pressure at the outlet is assumed constant and equal to zero.
One could claim that the application of the HT profile, used in CONF-1 and CONF-2, as an exit profile is an inappropriate choice since in many experimental investigations concerning round jets, see for example the work of Kyle and Sreenivasan (1993) or Boguslawski et al. (2013), the Blasius profile was found at the exit of the converging nozzle. However, as shown by Strykowski and Niccum (1991), the hyperbolic profile matches well the velocity distribution formed close to the nozzle exit (see Figure 4 in Strykowski and Niccum 1991). In addition, the numerical analysis of the global oscillations in low density jets Foysi et al. 2010;Boguslawski et al. 2016) demonstrated that the use of the HT profile gives reliable results in comparison to the experiments. Hence, it seems that this profile is only a slight simplification which, on the other hand, allow us to compare the LES predictions with the results of the linear stability theory directly. The impact of the velocity profile was discussed for the configuration CONF-2 in Sect. 4.2 for the completeness of the analysis.

Flow and Computational Parameters for CONF-1 and CONF-2
The Reynolds number in the simulations of the flows in CONF-1 and CONF-2 defined as Re D 1∕2 = U 1 D 1∕2 ∕ , where the symbol denotes the kinematic viscosity, was equal to Re D 1∕2 = 5000 . The analysis was performed for various velocity ratios I = 0.0 − 0.3 and shear layer thicknesses characterised by the parameter D 1∕2 ∕ = 24, 32, 40 . Choosing a relatively thick shear layer avoids interactions with the self-sustained convective oscillations reported by Boguslawski et al. (2013Boguslawski et al. ( , 2019 which appear for a shear layer characterized by D 1∕2 ∕ > 50 . It should be stressed that such thick shear layers as those considered in the paper are not typical for high Reynolds numbers but possible to accomplish by applying cylindrical nozzle tips (Pawlowska and Boguslawski 2020). The ones used in (Pawlowska and Boguslawski 2020) had a length of 1 − 25D and allowed the measurements of jets with the Reynolds numbers in the range from 5000 up to 20,000 and D∕ = 24 − 162 . The governing parameters of the test cases analysed for CONF-1 and CONF-2 are gathered in Table 2. Particular cases are denoted as C1/C2-D 1∕2 ∕ -I, where Mean inlet velocity profiles in CONF-3 assumed at the inlet plane of the main nozzle and in the exit section of the opposite nozzle C1 and C2 stand for CONF-1 and CONF-2, respectively. The computations were carried out on the mesh consisting of 272 × 110 × 64 points in axial, radial and azimuthal directions, respectively. Preliminary computations showed that such a mesh density is sufficient to obtain virtually grid independent results. Figure 4 shows the locations of the grid points across the jet region and in the shear layer for the velocity profiles with D 1∕2 ∕ = 40 . Assuming a shear layer thickness ( 99 ) as the region where 0.01U 2 < U m (r) ≤ 0.99U 1 the number of nodes across 99 depends on D 1∕2 ∕ and for D 1∕2 ∕ = 24, 32, 40 it is equal to 16, 12 and 10, respectively.
One of the necessary conditions for the occurrence of the global instability is a low level of the inlet turbulence intensity T i = u turb ∕U 1 , as shown by Monkewitz et al. (1990). The existence of too high turbulence intensity tends to inhibit the global instability. In all the simulations we assume T i = 0.1% , which corresponds well to the values of the turbulence intensity reported in the experimental researches of Strykowski and Wilcoxon (1993), where T i was less than 0.15% in the jet centreline and achieved maximum ( 1.5% ) in the shear layer; Asendrych (2007) with T i < 0.5% across the whole jet region; Favre-Marinet and  and Boguslawski et al. (1999), where T i = 0.2% was observed in the jet centreline and increased only in the shear layer region to T i = 8% . Regarding differences between the present settings and the experimental values of T i an important observation was made by Strykowski and Wilcoxon (1993) who showed that above I cr the level of T i did not have a substantial influence on the jet development.

Flow and Computational Parameters for CONF-3
The design of CONF-3 was aimed at the elimination of the suction of the fluid from the main jet and the formation of the counter-current mixing layer with a longer section of the parallel base flow. Figure 6a presents the iso-surfaces of instantaneous values of the Q-parameter with the axial velocity contours in the main cross-section. The Q-parameter is commonly used to indicate organised vortical structures and is defined as is the vorticity tensor. The embedded figure shows the velocity vectors and the vorticity contours in the region between the nozzles. These results correspond to the initial phase of the computations when a quasi-laminar base flow is formed between the nozzles and the transition to globally unstable regimes have not happened yet. The vortical structures appear only downstream, inside the opposite nozzle. As will be shown later when the global instability phenomenon begins the solution looks significantly different. For the present moment, it can be seen that CONF-3 meets the assumed goal. Especially, if one compares CONF-3 and CONF-1 for which the velocity contours presented in Fig. 2 are deflected right after the nozzle exit. In CONF-3 the region of the parallel flow is well formed, it extends over a longer distance and, what is most important, there is no suction of the fluid from the main jet stream.
The Reynolds number in CONF-3 is defined as Re D = U IP D∕ and equals to 5000. The Blasius velocity profile assumed at the inlet plane changes when the flow evolves inside the nozzle. The near-wall boundary layer thickness and the centreline velocity slightly grow. When the jet issues from the nozzle the velocity profile develops further and, as can be seen in Fig. 6b, in a short distance it transforms to an almost perfect HT profile applied in CONF-1 and CONF-2. Note that the locations x∕D = 1.2 and x∕D = 1.6 at which the profiles are presented correspond to the distances x∕D = 0.2 and x∕D = 0.6 behind the main nozzles.

3
The jets in CONF-3 were analysed only for two D 1∕2 ∕ values, D 1∕2 ∕ = 32 and D 1∕2 ∕ = 40 , and for I in the range of I = 0.0 − 0.225 . Increasing the I parameter further caused an amplification of a reverse flow in the upper nozzle, narrowing the shear layer of the main jet and strong asymmetry of the annular jet. To match the required D 1∕2 ∕ and I the values of U AI at the inlet to the upper nozzle (see Fig. 3) and D∕ defining the Blasius profile at the inlet to the main nozzle were adjusted by a trial-and-error approach. This led to D∕ = 48 and D∕ = 72 and U AI presented in table 1. The symbols in Fig. 6b represent the LES results at the time moment prior to the transition after which the flow becomes strongly unsteady and turbulent. It can be seen that the symbols correspond very well to the assumed HT profile. A small over-prediction of the velocity at r∕D 1∕2 > 0.8 is caused by the imposed co-flow velocity.
After the adjustment of U AI and D∕ the centreline velocity at the exit from the main nozzle was equal to U 1 = 1.02 − 1.04U IP and D 1∕2 = 0.91 − 0.94D that resulted in

Linear Stability Analysis
In this section, we apply the spatio-temporal analysis of the absolute instability phenomenon that will allow a proper interpretation of the LES results presented in the remaining part of the paper. The stability calculations are performed for the shear layer thicknesses corresponding to these used in the LES predictions. Although a parallel flow model used in LST is a significant simplification compared to the real flow conditions, the fundamental outcomes of this theory were confirmed by experimental studies of Strykowski and  Niccum (1991) and Strykowski and Wilcoxon (1993). We will interpret the LES results through the LST solutions as proposed by Jendoubi and Strykowski (1994). The small disturbance of flow parameters is assumed as a wave travelling along axial x and azimuthal directions with the amplitude varying along the radial direction r in the following form: where k is the complex wavenumber, is the complex frequency and m is the azimuthal mode number ( m = 0 represents the axisymmetric mode considered in the paper). By introducing equation 7, as shown by Jendoubi and Strykowski (1994) the equations of motions of small disturbances can be combined together to yield a second order ordinary differential equation in p(r) , thus: where c = ∕k is the complex phase velocity, Ma is the Mach number and S stands for density ratio. In the present paper only incompressible ( Ma = 0 ), homogeneous density ( S = 1 ) jets are considered. Hence, the stability equations for p(r) is simplified to: In the jet region where dŪ∕dr tends to zero, which is the case close to the jet axis and very far from it, Eq. 10 has the following asymptotic solution: where C 1 and C 2 are arbitrary constants which must be determined from the boundary conditions, and I m and K m are the modified Bessel functions of order m. Taking into account the following properties of the modified Bessel functions: and keeping in mind that a pressure perturbation amplitude is finite at the jet axis and tends towards zero at the jet periphery, the following asymptotic boundary conditions can be formulated: The stability equation was solved by the spectral tau method with an approximation of the eigenfunction of the pressure perturbation by the series of Chebyshev polynomials Canuto et al. (1988). The base flow is defined in terms of the mean velocity Ū (r) through the HT profile defined in Eq. (3) without the co-flow, i.e., The analysis is performed for D 1∕2 ∕ varying from 20 to 100 and R corespondig to the velocity ratio I = 0 − 0.3 The obtained solutions are in a form of complete maps of (k) ( -complex frequency, k-complex wave number) for which the branch point 0 is precisely found using an iterative algorithm proposed in Monkewitz and Sohn (1988) involving the shooting method.

Absolutely Unstable Mode I and Mode II
In agreement with the results of Jendoubi and Strykowski (1994) the stability analysis reveals two absolutely unstable modes. Figure 7 shows the absolute-convective instability boundaries for Mode I and Mode II as a function of D 1∕2 ∕ . It can be seen that in the entire range Mode II appears for lower values of I, for instance for D 1∕2 ∕ = 40 Mode II occurs at I cr = 0.088 while Mode I at I cr = 0.158 . The evolution of I cr with D 1∕2 ∕ for particular mode are significantly different. In the case of Mode I, initially for 20 < D 1∕2 ∕ < 24 I cr increases, it takes the maximum I cr = 0.184 at D 1∕2 ∕ = 24 and then smoothly drops down. In the case of Mode II, initially in the range D 1∕2 ∕ < 40 the values of I cr decrease reaching the minimum I cr = 0.088 at D 1∕2 ∕ = 40 . Then I cr smoothly rises towards the values of I cr found for Mode II. As shown in Jendoubi and Strykowski (1994) for D 1∕2 ∕ >> 100 the values of I cr are nearly the same for both modes and approximately equal to 0.14. Figure 8 presents the Strouhal number of Mode I and Mode II defined as

Characteristic Frequencies and Perturbation Wave Lengths of Mode I and Mode II
It can be seen that St D 1∕2 of both modes grow with I and D 1∕2 ∕ , but in the case of Mode I for larger D 1∕2 ∕ it reaches significantly higher levels. St D 1∕2 of Mode I is rather weakly dependent on I but it visibly rises with increasing D 1∕2 ∕ . Similar dependence of St D 1∕2 is found also for Mode II. In this case the impact of D 1∕2 ∕ is also small, it causes only a slight growth of St D 1∕2 , especially for D 1∕2 ∕ > 60 . It should be stressed that in this range the St D 1∕2 of Mode II is roughly two times lower than in the case of Mode I. Figure 9 presents the perturbation wave length ∕D 1∕2 of Mode I (a) and II (b), where = D 1∕2 ∕Real(k 0 ) . Unlike the St D 1∕2 profiles, which reveal similar behaviour for Mode I and Mode II, the evolution of ∕D 1∕2 shows the opposite trends. As seen in Fig. 9 ∕D 1∕2 of Mode I decreases with D 1∕2 ∕ and increases in the case of Mode II. For Mode I the profiles of ∕D 1∕2 overlap when D 1∕2 ∕ = 20 − 30 for all I values, but for higher values of D 1∕2 ∕ the influence of I is pronounced, the higher I is the shorter is the wavelength. In the case of Mode II the impact of I is more or less constant over the entire D 1∕2 ∕ range. It can be seen that ∕D 1∕2 decreases with the increase of I.
In general, it can be concluded that I weakly affects St D 1∕2 and ∕D 1∕2 whereas changes of D 1∕2 ∕ have large impact on these parameters, especially for D 1∕2 ∕ < 60.

Spatial Distribution of the Pressure Perturbation
A clear difference between two modes is seen in Fig. 10 showing the radial distribution of the magnitude of the pressure perturbation. Its values are normalised by the maximum values as |p|∕|p max | . As can be seen in Fig. 10b in the case of Mode II the perturbation extends across the whole jet for all I and D 1∕2 ∕ values. Except the results for D 1∕2 ∕ = 24 its dependence on I is small and in all but one cases the maximum of |p|∕|p max | occurs at the centreline. For D 1∕2 ∕ = 24 the increase of I causes the amplification of the perturbations in the shear layer close to r∕D 1∕2 = 0.5 . For I = 0.3 the maximum of |p|∕|p max | rises even above its value at the centreline. In some sense the behaviour of Mode II for D 1∕2 ∕ = 24 and I = 0.3 is similar to that observed for Mode I. As can be seen in Fig. 10a in this case the trend of |p|∕|p max | exhibit a completely different behaviour. The profiles obtained for D 1∕2 ∕ = 24 have the global maximum at the centreline and the local maximum in the region of the shear layer. Increasing D 1∕2 ∕ to 40 or 100 causes lowering of |p|∕|p max | at the centreline and its sudden jump near r∕D 1∕2 = 0.5 . For D 1∕2 ∕ = 40 the impact of I can be regarded as significant and larger than for D 1∕2 ∕ = 100 for which the maxima of |p|∕|p max | are concentrated around the shear layer and vanish at the jet periphery as well as at the jets axis.
Summing up, the amplitude of the perturbation related to Mode I visibly depends on D 1∕2 ∕ and for its smaller values also on I. In the case of Mode II the impact of D 1∕2 ∕ is definitively smaller whereas I affects the perturbation only when D 1∕2 ∕ is small.

Results of LES
The results obtained from the spatio-temporal stability analysis will be used as a qualitative and quantitative indicator allowing us to distinguish the stability mode predicted by LES and to assess to what extent the obtained solutions are consistent with the theoretical basis. Note that the perfect agreement cannot be expected as the LES are performed for real viscous fluids while the results of stability analysis are obtained for an inviscid flow. In all the cases the scenario of the simulations was the same. The jets were injected impulsively and then tracked as it developed in the domain. The calculations were conducted with the constant time step t = 0.00125D 1∕2 ∕U 1 . After the time 100D 1∕2 ∕U 1 , when the flow was fully developed, the time-averaging procedure started for the next 300D 1∕2 ∕U 1 which corresponded to 20 flow through times. In these cases the fluctuations start growing at x∕D 1∕2 = 2.0 and attain the maximum level of 20% of U 1 at the distance approximately equal to 9D 1∕2 . A faster perturbation growth is observed starting from x∕D 1∕2 ≈ 1 for I = 0.15 − 0.275 . Additionally, for I = 0.25 and 0.275 the profiles of fluctuations are characterised by the local maxima that correspond to the local minima of the mean velocity profiles. A sudden change is found for I = 0.3 . In this case, the fluctuations grow right from the inlet and attain two maxima. The first one on the level of 10% of U 1 at the axial distance x∕D 1∕2 = 2.0 and the second one further downstream where the flow is fully developed. It should be stressed that the second maximum is located closer to the inlet and achieves a lower value than in the other cases with lower I.  above which the absolute mode emerges. The increase of D 1∕2 ∕ to 40 changes I cr which in this case is equal to 0.275 (see Fig. 11c). Moreover, one should notice that the first local maxima for the cases with I ≥ I cr are lower than in the cases with D 1∕2 ∕ = 24, 32 . The fluctuations attain the level of 5% and 4% of U 1 for I = 0.275 and 0.3, respectively. Similar damping of the velocity fluctuations was observed by Strykowski and Niccum (1991). In their studies this effect was even stronger, however, they analysed a jet with a considerably thinner shear layer D∕ = 230 . Taking into account the impact of I on the evolution of the perturbation along the jet axis discussed in Sect. 3.3 it can be readily concluded that the observed flow behaviour is typical for Mode I for which the increase of I reduced the perturbation at the centreline, and this effect was stronger for larger D∕ . Figure 12 shows the profiles of the fluctuations in a logarithmic scale along the x-axis at the radial distance r∕D 1∕2 = 0.5 , which corresponds to the shear layer region, where according to the spatio-temporal stability analysis the clear manifestation of a change in  Fig. 13. As can be seen, the differences between these two regimes are very clear. Contrary to the intermediate regime, where the small vortex forms at the axial distance ≈ 1.4D 1∕2 , in the global regime two stronger vortical structures are visible, the first one close to x∕D 1∕2 ≈ 0.5 where the profile of fluctuations attained the maximum value (see Fig. 11c). Downstream, the vortices break up within a distance smaller than x∕D 1∕2 = 2.0. frequency for the Kelvin-Helmholtz regime is presented in Fig. 14a. One can observe that in these cases the spectra are characterised by a broadband distribution centred at St D 1∕2 ≈ 0.6 for C1-24-0.150 and St D 1∕2 ≈ 0.7 for C1-32-0.150 and C1-40-0.150. A similar distribution is found for the intermediate regime, however, in this case the broadband peaks are present already close to the nozzle at x∕D 1∕2 = 1.0 . In the global regime the spectra presented in Fig. 14c reveal existence of precisely defined distinct peaks which confirms the triggering of the globally unstable mode. The peaks occur at St D 1∕2 = 0.75 , 0.91, 0.96 and 1.02 for the cases C1-24-0.300, C1-32-0.300, C1-40-0.275 and C1-40-0.300, respectively. In the comparison with the linear stability analysis (see Fig. 8) the characteristic St D 1∕2 are higher, however, the most important observation is that exactly the same trend is observed. The peaks move towards larger St D 1∕2 with the increasing D 1∕2 ∕ and I. It should be mentioned here that even in DNS results of  of low density jets the global mode frequency was significantly higher than the absolute mode frequency obtained from the linear stability theory.

Summary of the Simulations for CONF-1
Despite the fact that the LES and linear stability results are in good qualitative agreement, one should keep in mind that the critical velocity ratios were found to be significantly larger. The linear stability shows I cr considerably lower than 0.2 while in the LES they are I cr = 0.3 for D 1∕2 ∕ = 24 and D 1∕2 ∕ = 32 and I cr = 0.275 for D 1∕2 ∕ = 40 . The critical velocity ratio higher than the theoretically predicted was also observed by Strykowski and Niccum (1991), who analysed the influence of the shear layer thickness on the global mode for the jets with D∕ in a range 135 − 230 . They showed that the decrease of D∕ from 230 to 135 increases I cr from 0.14 to 0.18, approximately. Having in mind that the present analysis is performed for significantly thicker shear layers, for which LST shows the fast growth of I cr when the D 1∕2 ∕ becomes lower, the obtained results from LES can be certainly treated as accurate and showing the existence of global Mode I.

LES of the Counter-Current Jets for CONF-2
An analogous analysis was conducted for the second configuration (CONF-2) for which a tube was used to extend the region of the parallel counter-flow. As seen in Fig. 15 where the white isolines denote the zero axial velocity, an application of the wall around the main jet causes a significant extension of the counter-current mixing layer and prevents hitting the jet by the surrounding flow. In the case of CONF-2 the results obtained are compared with the results of spatio-temporal stability analysis and also with the experimental data of Favre-Marinet and , who studied the jet with D∕ ≈ 40 . The solutions for these conditions are presented in Figs. 16 and 17. They were obtained applying the HT and Blasius inlet velocity profiles. It can be seen in Fig. 16 that for the jet with I = 0.0 the results obtained with the HT profile are in a reasonable agreement with the experimental measurements. Although the length of the potential core and initial growth of the fluctuations at x∕D 1∕2 < 6 agree quite well it can be seen that the occurrence of the maximum of fluctuations is predicted approximately 2.5D 1∕2 further from the inlet and its level is nearly 30% higher compared to the experiment. Comparing the profiles predicted for jets with I = 0.2 and 0.3 with the experimental data one can observe that the most visible differences appear for I = 0.3 . The decay of the mean velocity is faster and the fluctuating velocity profiles attain a higher maximum level. Nevertheless, the observed discrepancies have only quantitative character. For I = 0.2 the profiles match almost perfectly. Although the HT profile is commonly applied in the jet flow simulations, it should be stressed that it corresponds rather to the profile formed downstream of the nozzle exit than to the profile at the nozzle exit. As presented by Wawrzak et al. (2015) the application of the Blasius profile, which better reflects the experimental velocity profiles, might give results closer to experimental data. Indeed, as can be seen in Fig. 17 the use of the Blasius profile improves the solutions, particularly for the cases with I = 0.2 and 0.3 whereas for I = 0.0 the profiles practically do not differ. The maximum of the fluctuations remains over-predicted compared to the measured values. Most likely this discrepancy is caused by different inlet turbulence characteristics (turbulent length and time scales), the influence of which on the jet behaviour in the near-field region is of crucial importance, particularly at low Ti level (Wawrzak et al. 2015). Contrary to the solution obtained for I = 0.0 the results obtained for I = 0.2, 0.3 turn out to be in excellent agreement with the experimental data. The mean velocity values are slightly under-predicted only for I = 0.3 , whereas the locations and levels of the fluctuations maxima are captured very accurately in both cases. We attribute this to the fact that above the critical velocity ratio the initial turbulence characteristics have only a minor impact on the jet development (Strykowski and Wilcoxon 1993). Despite the fact that the application of the Blasius profile leads to more accurate results confirming the correctness of the performed LES in the following part of the paper we stick to the HT profile as it enables a direct comparison with the results of the spatio-temporal analysis.
The mean and fluctuating velocity profiles along the jet axis for D 1∕2 ∕ = 24 and D 1∕2 ∕ = 32 that complements the results obtained for D 1∕2 ∕ = 40 in Fig. 16 are presented in Fig. 18. It can readily be seen that the profiles qualitatively differ from those obtained for CONF-1 in three main aspects: (i) a sudden change of the flow behaviour is observed already for I cr = 0.125 for D 1∕2 ∕ = 24 and I cr = 0.1 for D 1∕2 ∕ = 32 and D 1∕2 ∕ = 40 ; (ii) the mean profiles are characterised by the faster decay with an increase of I; (iii) the perturbation growth is observed right from the inlet plane and the maximum of the fluctuation increases with I. Taking into account the results of LST one could conclude that CONF-2 provides conditions favourable for triggering Mode II. As presented in Fig. 7, I cr of absolute Mode II predicted by the spatio-temporal analysis agrees quite well with I cr obtained in the LES. Also the trend of I cr in the function of D 1∕2 ∕ is accurately captured. Additionally, the LES results show that the impact of D 1∕2 ∕ on the solutions for I above I cr is negligible. Only for I = 0.3 a slightly faster decay of the velocity is noticeable for higher D 1∕2 ∕ values. Consequently, the maximum of the fluctuations is shifted upstream. Figure 19 shows the amplitude spectra of the axial velocity fluctuations registered at the jet axis and the distance x∕D 1∕2 = 1.5 for the cases with the velocity ratio I ≥ I cr . Contrary to CONF-1, for which the distinct peaks were observed, the present results reveal rather broadband distribution. In some cases, however, the characteristic frequency, which indicates the periodic formation of strong vortical structures can easily be found. In the case of  0.15, respectively. In general, similarly as for CONF-1, it was found that an increase of I or D 1∕2 ∕ leads to the higher St D 1∕2 values that is consistent with the linear stability analysis. It is worth noting that the enhancement of the fluctuations in a large frequency bandwidth without distinct peaks for larger I values was observed in Favre-Marinet and . It could be observed in Fig. 19 that for the thinner shear layer, D 1∕2 ∕ = 40 , the peaks in spectra are less distinct than for the thicker ones. This might seem somewhat surprising as the opposite tendency could be expected taking into account that the absolute growth rate increases for the thinner shear layer. Such a behaviour can be attributed to the fact that the perturbation, which wavelength is growing for thinner shear layer in the case of the mode II, cannot evolves within a limited range of parallel flow in CONF-2. Hence, as longer perturbation waves have no sufficient distance to develop their occurrence is less favourable that results in irregular flow oscillations leading to broadband peaks in the velocity spectra.

Velocity Spectra
We found that an increase of I (e.g. to I = 0.2 for D 1∕2 ∕ = 40 ) seems to completely destroy the peaks. Such a behaviour is in contradiction with the results obtained by Strykowski and Wilcoxon (1993) who observed peaks for I in a similar range. However, they analysed the jets with a significantly higher Reynolds number 34,000 ≤ Re D ≤ 110,000 and thinner shear layer D∕ > 100 . Moreover, Strykowski and Wilcoxon (1993) used an

3D Analysis of an Instantaneous Flow
The evolution of the spectra with the increase of I suggests that there must be a significant change in the flow regimes between the cases with lower and higher values of I. The differences are readily seen in Fig. 20 showing the iso-surfaces of the Q-parameter and the axial velocity contours for the cases C2-40-0.100 and C2-40-0.200. The  Fig. 20 show: (i) the axial velocity contours in a cross-section perpendicular to the jet axis at the distance 3D from the inlet; (ii) the vorticity contours and the velocity vector field in a vicinity of the inlet -the region marked with the black rectangle placed between the jet axis and the extension tube. The instantaneous values of Q-parameter reveal the formation of large vortical structures only for I = 0.1 (Fig. 20a). As can be seen in the vector distribution in this case the conditions for establishing almost perfectly parallel flow between the jet and counter-current stream are favourable. It can be seen that toroidal vortices form already very close to the nozzle exit ( x∕D 1∕2 ≈ 1.0 ). They successively grow up to the distance x∕D 1∕2 = 3.0 , mutually interact and eventually break up downstream. The axial velocity contours at the main cross-section reveal the occurrence of a separate stream inclined approximately 45 deg to the jet axis. This stream is released from the main jet core, it is called side-jet and was also observed in the experiments Monkewitz et al. 1990;Hallberg and Strykowski 2006). In the cross-section perpendicular to the jet axis it can be seen that there are five side-jets, which are non-uniformly distributed in the azimuthal direction. The mechanism of their generation is related to the interaction between the toroidal vortical structures. The side-jets are a manifestation of the global instability phenomenon, yet they cannot be treated as direct evidence of its occurrence. Regarding the impact of I Fig. 20b shows that the increase of the I to 0.2 completely alters the flow picture. The toroidal vortices do not exist anymore and the flow field seems to be highly turbulent. The vector field shows that the fluid from the main jet is taken to the suction slot and consequently the parallel region, which existed for I = 0.1 , is completely destroyed in this case. As it was already mentioned such a destructive influence of extensive suction that prevents the formation of strong vortices was observed in Boguslawski et al. (1999). One can conclude that the high velocity ratio suppresses or at least significantly disturbs the global instability phenomenon.

Summary of the Simulations for CONF-1 and CONF-2
It has been shown that the application of the suction around the main jet allows emerging of the global Mode I and Mode II. The differences between the configurations CONF-1 and CONF-2 were found to be mainly in the streamwise distances over which the region of parallel counter-flow extends. The length of this region undoubtedly impacts on the flow response. It is known that triggering of the global modes is possible when the counterflow is parallel at the distance of the order of the perturbation wave length (Jendoubi and Strykowski 1994). As could be seen in Fig. 8 Fig. 15 the region of the counter-flow for CONF-1 exists only up to the axial distance x = 0.3D 1∕2 , while for CONF-2 it extends to x = 0.65D 1∕2 . This means that in CONF-1 the likelihood of the occurrence of Mode I is definitively larger, whereas the conditions favourable for Mode II are maintained in CONF-2 where the region of the parallel oriented flow extends further downstream. Consistently with the above expectations, CONF-1 allows triggering the Mode I and CONF-2 facilitates the appearance of Mode II. However, it should be stressed that a significant increase of I in CONF-2 leads to the suppression of the global mode due to sucking the fluid from the main jet stream. To overcome this problem the configuration CONF-3 is introduced in which this effect does not have place and the parallel flow region extends over a longer distance.

LES of the Counter-Current Jets for CONF-3
Configuration CONF-3 presented in Fig. 3  vertical line at x∕D = 1.0 indicates the end of the main nozzle, which in CONF-3 is the part of the computational domain. Regarding the impact of I it can be seen that below I = 0.15 the profiles are only slightly different from the case with I = 0.0 . For I < 0.15 the fluctuations start growing at the distance x∕D ≈ 2.0 (i.e. 1.0D from the nozzle exit) and attain the level around 0.2U 1 at x∕D = 6.0 − 8.0 . The situation changes drastically for I ≥ 0.15 as in these cases two local maxima are readily seen. The first one is located at the distance x∕D ≈ 2.5 and it reaches the level 0.06U 1 − 0.12U 1 depending on I and D 1∕2 ∕ . Hence, one can assume that I = 0.15 is the critical value above which the global instability emerges. The occurrence of the fluctuation maximum close to the inlet is characteristic of a rapid formation of strong vortices, that in the present case exist only over a short distance. They break down flowing inside the opposite upper nozzle that is accompanied by a sudden drop in the fluctuations in between x∕D = 2.5 − 4.0 . Further downstream the fluctuations intensities recover and grow up to the level 0.18U 1 . The second maxima occur in more or less the same locations as for the cases with I < 0.15 . The configuration CONF-3 is certainly much more difficult for analysis and interpretation due to the geometry complexity. One could suspect that the presence and the close proximity of the two nozzles lead to the interaction of each of the jets with the opposite nozzle. However, a close examination of the flow field does not show any stagnation zones indicating such an interaction. The walls of the opposite nozzle only impose some constraints on the development of vortical structures travelling downstream, however, their interactions with the walls do not impact the mechanism of their generation. Excluding an influence of the walls on the instability scenario one could, however, still suspect that the dynamics of two opposite jets is determined by an interaction of convective waves travelling in the opposite directions. It should be stressed that analysis of the flow development showed that the flow transition starts right after the laminar steady base flow is established with a sufficiently strong counter-current velocity. At this stage no convective waves were observed. Hence, despite the geometry complexity we are convinced that interpretation of the LES in terms of the linear theory like in CONF-1 and CONF-2 is justified. Trying to identify whether the observed instability is characteristic for Mode I or Mode II we analyse dependence of the results on D 1∕2 ∕ and I both qualitatively and quantitatively. In Fig. 21 it can be seen that: • for I > I cr = 0.15 the maxima of the fluctuations decrease with the increase of D 1∕2 ∕ .
According to LST (see Fig. 7a) such a behaviour is typical for Mode I and it was observed also in CONF-1 (see Fig. 11) in which this instability mode was found. • for I > I cr = 0.15 the maxima of the fluctuations increase with the increase of I. Such a behaviour was found in CONF-2 (see Fig. 18) characterised by Mode II, though in CONF-2 the impact of I was much more pronounced. • for I > I cr the local maxima occur close to the nozzle exit and their locations and levels are very close to that found in CONF-1. • I cr = 0.15 obtained for CONF-3 is close to I cr ≈ 0.165 of the global Mode I in the range D 1∕2 ∕ = 30 − 40 , as predicted by LST (see Fig. 7), and also close to I cr = 0.14 reported by Strykowski and Niccum (1991).
Taking into account the above findings it can be said that the only result indicating the occurrence of Mode II in CONF-3 is the qualitative similarity of the dependence of the fluctuations maxima on I. Hence, we are inclined to state that the global instability observed in CONF-3 is characteristic for Mode I. Additional convincing arguments supporting this supposition are provided in the following subsection.

Velocity Spectra
The amplitude spectra of the axial velocity at the location x∕D = 1.5 (0.5D above the main nozzle) are shown in Fig. 22. Similarly as in CONF-1 (see Fig. 14 in LES agree very well with that found in LST, while for the remaining cases they are larger, especially for D 1∕2 ∕ = 32 . An over-prediction of the theoretical St D 1∕2 values was previously reported in  and Boguslawski et al. (2016) in DNS and LES studies of low density jets. It seems that such a behaviour is typical in numerical simulations and can be attributed to non-linear interactions of large flow structures, not present in LST analysis, that cause the increase of St D 1∕2 . Beside this discrepancy the influence of D 1∕2 ∕ and I on the variation of St D 1∕2 is the same in LES and LST. It can be seen that St D 1∕2 rises visibly with the increasing D 1∕2 ∕ and only slightly with I.

3D Analysis of an Instantaneous Flow
The occurrence of distinct peaks in the spectra in between the nozzles, i.e. at x∕D = 1.5 , suggests that the origin of the global instability is located in the counter-current mixing region. Figure 24 shows the iso-surfaces of the Q-parameter and the axial velocity contours at two different time instances for the case with I = 0.15 and D 1∕2 ∕ = 32 . The first moment presented in Fig. 24a corresponds to an initial phase of the flow development when the toroidal vortices are formed far from the inlet and when the mixing layer between the jet and the counter-current flow is nearly axisymmetric and seems to be laminar. The counter-current flowing streams are almost perfectly parallel that implies suitable conditions for the development of the absolute instability. Indeed, this is what happens shortly after. A strong vortical structure forms between the two nozzles, which in Fig. 24b exhibits as the toroidal vortex associated with the vorticity maximum at ( r∕D = 0.45 , x∕D = 1.7 ). In the later time, the vortices in this region appear continuously thus causing the self-sustained oscillations leading to the global instability.

Conclusion
The paper presented the analysis of the global instability phenomenon in the counter-current round jets. The research was performed using the LES method, the accuracy of which was confirmed by the comparison with experimental data, and referring to previous works where the numerical code applied was used and turned out to be accurate. We studied two configurations (CONF-1, CONF-2) of the round jet with the suction applied around the nozzle to establish the counter-current flow in the same way as in earlier experimental researches (e.g. Strykowski and Niccum 1991;Strykowski and Wilcoxon 1993;Jendoubi and Strykowski 1994;Boguslawski et al. 1999). Additionally to these typical set-ups, we proposed, a new configuration (CONF-3) in which the counter-current flow resulted from the fluid streams released from two opposite nozzles. The control parameters examined in the research were the velocity ratio I and the shear layer thickness characterised by the D 1∕2 ∕ parameter. The presented results univocally revealed the occurrence of the global instability phenomenon in all the configurations analysed. It was shown that the type of the instability mode (Mode I and Mode II) and the range of I at which it appeared were dependent on the configuration applied. The simulations performed for CONF-1 revealed that in this case the jets exhibit behaviour typical for Mode I characterised by two maxima in the profiles of the axial velocity fluctuations and sharp peaks in the spectra of the axial velocity. The range of I for which the flow manifested the global instability regime was very narrow and weakly depended on the shear layer thickness. The region of the parallel flow required to trigger the global instability phenomenon was very short and existed only for I cr larger than predicted by LST. Attempts of increasing I cr resulted in the suction of the fluid from the main jet and destruction of the parallel mixing layer caused by the fluid flowing radially towards the jet from the ambient region. The fact that the counter-current streams were parallel only over a short distance excluded the possibility of occurrence of Mode II for which this distance had to be definitively longer, as shown by LST.
The aim of the extension tube used in CONF-2 was to eliminate the radial flow from the surroundings and to form a longer section of parallel counter-current streams at lower I. The results obtained showed that CONF-2 fulfils these goals. The value of I cr causing the emergence of the global instability regimes was in good agreement with LST and, as expected, the Mode II appeared. However, the range of I > I cr for which this happened was limited as the extension tube did not prevent the sucking of the fluid from the main jet, which, similarly as in CONF-1 had a destructive influence on the global instability phenomenon.
In CONF-3 none of the negative effects mentioned above took place. The external stream issuing from the upper nozzle smoothly flew along with the main jet and when it approached the lower nozzle it was oriented radially in the direction outside due to the inclined shape of the nozzle. In such conditions the mixing layer was almost undisturbed and the counter-current flow region was parallel over the distance comparable to that observed in CONF-1. Such a flow picture was, however, only temporary. Small disturbances developed when the simulation time passed and at a certain time moment the flow became strongly unsteady and turbulent, and for I ≥ I cr switched to the globally unstable regime. In this case, both the value of I cr and its dependence on D 1∕2 ∕ , as well as the locations of the peaks in the velocity spectra in the function of D 1∕2 ∕ agreed well with those for Mode I obtained from LST and in previous experimental investigations. The conditions facilitating the occurrence the global instability in CONF-3 were found to be also limited to some range of I, which substantial increase caused asymmetry and turbulization of the flow inside the gap of the upper nozzle leading to the destruction of the counter-current parallel flow region.
In summary, trying to respond to the questions formulated in the introduction it can be said that in counter-current jets: (i) for I significantly larger than I cr the global instability disappears due to the suction of the fluid from the main jet and lack of the parallel flow region between the counter-current streams, as in CONF-1 and CONF-2, or because of high turbulence intensity in the external stream, as in CONF-3; taking this into account it seems that the control of the global instability phenomenon in a wide range of ≥ I cr is hardly possible; (ii) having the flow conditions favourable for triggering the global instability, i.e., a low level or no suction of the fluid from the main jet and the parallel counter-current base flow, the type of the appearing instability mode (Mode I / Mode II) depends primarily on the length of parallel flow region; when this is short only Mode I has the chance to emerge, as in CONF-1 and CONF-3. Having in mind the geometry of CONF-3 it seems possible that by changing the axial distance between the nozzles the length of the parallel flow could be controlled allowing the occurrence of Mode I when the nozzles are close to each other and I is large (I > I cr,Mode I ≈ 0.15) or Mode II when the nozzles are spaced farther from each other and I is around 0.1, i.e., I cr,Mode II < I < I cr,Mode I . Potentially, such a flow control technique can find applications in practical devices, e.g,. in novel MILD-type combustion chambers where an intense mixing of fresh reactants with recirculating combustion products in a counter-current stream is of crucial importance. This problem is left for future study.