4.1. General description

A typical flow picture displaying the liquid–gas interface and the velocity magnitude is presented in figure 1(b). Initially flat, the interface progressively deforms into a wavy shape. The liquid strongly interacts with the high-speed jet in the downstream part of the flow: the liquid jet gets thinner due to the transfer of momentum from gas to liquid and flow rate conservation, but one can also see large amplitude waves deviating the gas jet, as well as ligaments and liquid fragments pulled out by the jet. The gas–gas mixing layer shows the triggering of a Kelvin–Helmholtz instability. No recirculation is observed in the top part of the domain, showing that the coflow actually provides the mass entrained by the high-speed jet. The vertical position of the interface is sampled for all positions downstream of the inlet at a frequency of 10 000 Hz, allowing the study of the spatiotemporal development of the instability.

In figure 2(a), we show the spatial evolution of the wave amplitude, and in figure 2(b) the spectrogram of the interface height for case A2, with $n=2$ and $\delta _d=1$. The wave amplitude is computed using the method presented in Matas et al. (Reference Matas, Marty and Cartellier2011): for a given downstream position, we construct a histogram of the interface positions over the sampling time and exclude the lowest and highest $0.5\,\%$. The remaining width of the histogram is taken as the amplitude $A$ of the waves. Several zones can be distinguished. Focusing on figure 2(a), one can see that for $20 \leq x/\delta _g \leq 50$, the computed amplitude ($\circ$, blue) seems to increase exponentially with downstream distance,

(4.1)\begin{equation} A= A_0 \textrm{e}^{k_i x}, \end{equation}

with, for figure 2(a), $k_i=k_{i, exp}=323.2$ m$^{-1}$, where $k_{i, exp}$ is the value obtained experimentally (Matas et al. Reference Matas, Marty and Cartellier2011), as represented with a dashed line in figure 2(a).

Figure 2.

Spatiotemporal development of the instability for the case A2 with $n=2$, $\delta _d=1$. (a) Spatial evolution of the wave amplitude with: computed wave amplitude ($\circ$, blue); (4.1), with $A_0/\delta _g=3.8\times 10^{-4}$ and $k_i=k_{i,exp}=323.2$ m$^{-1}$ (orange dashed line); (4.2), with $C_0=0.21$ (orange solid line). The insert shows the same data in a log scale, with a close-up of the exponential growth region. (b) Spectrogram of the interface height. $P*$ is the normalized spectral power.

This observation has to be taken cautiously. The stability analysis predicts that the instability is absolute for these conditions, as seen on figure 3, where we see the pinching between the confinement branch, located close to the $k_i$ axis, and the shear branch. Nonlinearity is therefore expected to quickly dominate, as the instability will grow exponentially in time and space. As the group velocity at the pinching point is zero, the predicted spatial growth rate has no clear meaning for the waves travelling downstream of the injector (Huerre & Monkewitz Reference Huerre and Monkewitz1990).

Figure 3.

Spatial branches for case A2. (a) Pinching between confinement and shear branches, at $\omega _i=30$ s$^{-1}$. The confinement branch is located close to the $k_i$ axis. The arrows indicate the direction in which $\omega _r$ increases. (b) Zoomed-out view of the pinch point showing the behaviour of the shear branch at low and large $k_r$ for $\omega _i=30$ s$^{-1}$, and for two additional values of $\omega _i$. At large $\omega _i$, the shear branch is located in the $k_i>0$ half-plane, while the confinement branch remains close to the imaginary axis in the $k_i<0$ half-plane: the two colliding branches therefore originate from opposite half-planes, in agreement with the Briggs criterion.

However, as said in the introduction, the predicted spatial growth rate controls the scale of oscillations in the cross-stream direction. The confinement branch is crossing the $k_i$ axis around $-n{\rm \pi} /(2L)$, with $L$ the cross-stream characteristic length and $n=\{1,2,3,\ldots \}$. At the pinching point shown in figure 3, $k_{i,LSA} = - 180$ m$^{-1}$, which is close to $-{\rm \pi} /(2H_l)=-157$ m$^{-1}$, showing the matching between cross-stream perturbations and cross-stream confinement that leads to resonance, i.e. to absolute instability.

An exponential growth zone has been seen experimentally and numerically (Agbaglah et al. Reference Agbaglah, Chiodi and Desjardins2017; Ling et al. Reference Ling, Fuster, Tryggvason and Zaleski2019). Figure 2(a) shows that the wave growth obtained numerically is in agreement with the wave growth observed experimentally. The insert in figure 2(a) shows that the zone of exponential growth rate is extremely short, as seen in experiments.

The end of the previous zone occurs when the amplitude is of the order of the vorticity thickness, as observed in experiments. For $x/\delta _g \geq 75$, the amplitude is growing linearly with downstream distance, following well the self-similar model of Hoepffner, Blumenthal & Zaleski (Reference Hoepffner, Blumenthal and Zaleski2011),

(4.2)\begin{equation} A = C_0\sqrt{\rho_g/\rho_l}U_g \tau, \end{equation}

with the characteristic time taken as $\tau =x/U_D$, as in Marty (Reference Marty2015), $U_D$ is obtained using (1.2) and $C_0$ is a model constant to be specified. This model is represented by the continuous line in figure 2(a).

In figure 2(b), one can see the spatial evolution of the interface height spectra. For each location downstream of the inlet, we obtain the interface height spectrum and then normalize the spectral power by its local maximum. Note that we use the Welch method (Welch Reference Welch1967) to facilitate the extraction of the most amplified frequency, i.e. before taking its Fourier transform, we split the signal in 10 parts with a 60 % overlap between them, the resulting spectra are then averaged. A very clear frequency peak is observed and remains unchanged over a large extent of the domain, which may suggest the presence of a global mode of instability. Close to the outlet, lower frequency peaks are appearing.

In the next section, we detail the influence of the numerical parameters and the integration time on these results.

4.2. Assessment of the numerical methods

In figure 4(a), we show the influence of integration time on the most amplified frequency obtained at three downstream locations for two sets of injection conditions. The most amplified frequency is presented as a function of the normalized duration of the signal $\Delta T f_{th}$, with $\Delta T$ the signal duration and $f_{th}$ the predicted wave frequency obtained by linear stability analysis. Note that the statistics are collected once the averaged height is stabilized over time. One can see that as soon as $\Delta T f_{th}$ is larger than $10$, i.e. approximately 10 waves have been generated, the variations of the most amplified frequency in time and space are of the order of the spectral resolution. This is in agreement with the results of Agbaglah et al. (Reference Agbaglah, Chiodi and Desjardins2017). For both cases, the most amplified frequency is stable in time, and with the increase of spectral resolution all measured frequencies become perfectly independent of the position, providing that the latter is in the interval $25 \leq x/\delta _g \leq 150$. In order to reduce the uncertainty on the frequency value an integration time of $\Delta T f_{th}\geq 40$ would be ideal, but this choice will be more difficult to achieve for 3-D simulations, or even very resolved 2-D simulations. The present analysis shows that even an integration time of $\Delta T f_{th}\sim 10$ is sufficient to obtain a good estimate of the frequency.

Figure 4.

Dependence of frequency measurements on integration time and domain length for $n=2$, $\delta _d=1$. (a) Spatiotemporal convergence of the most amplified frequency with: Blue, case B1; magenta, case A2; $x/\delta _g=44$ ($\circ$); $x/\delta _g=67$ ($\triangle$); $x/\delta _g=89$ ($\square$). The error bars represent the spectral resolution. Horizontal dashed lines are the theoretical frequency for each case. (b) Dependence of frequency spectra on domain length for $\Delta Tf_{th}=40$, $x/\delta _g=89$ with: case B1 (Blue); case A2 (magenta); continuous line, $L_x/\delta _g=240$; dotted line, $L_x/\delta _g=320$.

In figure 4(b), we show the impact of domain length on the frequency spectra for two sets of injection conditions. Our tests show that, for the cases of figure 4(b), the most amplified frequency is not affected by the position of the outlet boundary condition once $L_x/\delta _g$ is larger than $240$, which has been confirmed for other cases (results not shown). We will therefore ensure that this condition is met for all cases.

In figure 5, we show the effect of resolution on the spatiotemporal development of the instability for case $B1$ and two values of the deficit parameter. In figure 5(a) is presented the evolution of the frequency spectra with $n$ ($=\Delta x/\delta _g=\Delta y/\delta _g$), the resolution, for $\delta _d=1$. We see that, as the mesh is refined, the most amplified frequency approaches the theoretical value. The error on the finest mesh, here $n=8$, is only $4.6\,\%$, while the error for the coarsest resolution, $n=2$, is $12\,\%$. In figure 5(b), the convergence of the amplitude with mesh resolution for $\delta _d=1$ is shown. From $n=2$ ($\triangle$, black) to $n=16$ ($+$, red) the amplitude is mostly affected by the resolution in the area close to the injector, where the amplitude is small. Far from the inlet, i.e. $x/\delta _g \geq 75$, the amplitude is essentially not affected by resolution and closely follows the self-similar model. Mesh resolution mostly affects the location at which the amplitude reaches this self-similar region. This location is shifted closer to the injector as the mesh is refined. Note that the amplitude is essentially not modified between $n=8$ and $n=16$, meaning that we reach independence from mesh resolution. Note also that, close to the injector, we do again see a zone of exponential growth, with a slope equal to the one observed experimentally for these injection conditions. As the mesh is refined, the transition between the inlet of the domain and the zone of exponential growth progressively disappears, which is not the case in experiments (Matas et al. Reference Matas, Marty and Cartellier2011) or simulations including a splitter plate (Agbaglah et al. Reference Agbaglah, Chiodi and Desjardins2017; Ling et al. Reference Ling, Fuster, Tryggvason and Zaleski2019).

Figure 5.

Effect of mesh resolution on the instability characteristics for case B1 and two values of the deficit parameter. In the zone of interest, $\Delta x = \Delta y = \delta _g /n$. All spectra are taken at $x/\delta _g= 67$, with $\Delta T f_{th} \geq 20$. (a) Convergence of frequency with mesh resolution for $\delta _d=1$. The vertical dashed line represents the corresponding theoretical frequency. (b) Convergence of wave amplitude with mesh resolution for $\delta _d=1$ with: (range dashed line), (4.1), with $A_0/\delta _g=0.016$ and $k_i=k_{i,exp}=569.9$ m$^{-1}$; (orange solid line) (4.2), with $C_0=0.21$. (c) Similar to (a), with $\delta _d=0.5$. (d) Similar to (b), with $\delta _d=0.5$, $A_0/\delta _g=0.046$ and $C_0=0.27$. Note that in (b,d), only one data point every six is shown to enhance readability.

In figure 5(c), we present the evolution of the frequency spectra with mesh resolution for $\delta _d=0.5$ and the same B1 case. Again, one can see that, as the mesh is refined, the most amplified frequency converges towards the predicted value, but the error on the finest resolution, $n=10$, is still $12\,\%$. The evolution of amplitude with resolution for $\delta _d=0.5$ is shown in figure 5(d). Similarly to the case without velocity deficit, we can see that as the mesh is refined, the position at which the amplitude enters the self-similar region is shifted towards the inlet. In the exponential growth region, we are again in very good agreement with the slope measured experimentally by Matas et al. (Reference Matas, Marty and Cartellier2011). Note that we do observe that the self-similar growth region is affected by the reduction of $\delta _d$ as we had to increase $C_0$ from $0.21$ to $0.27$ in order to obtain a good agreement between the computed amplitude and (4.2) in the downstream area. This is consistent with our observation of a decrease in the wave speed with the decrease of $\delta _d$ (results not shown). This suggests, as already pointed out by Marty (Reference Marty2015), that one needs to adjust the velocity scale used in (4.2) to the variations of the interfacial wave speed, therefore, to the details of the velocity profile, in order to obtain a good self-similarity of the wave amplitude.

Overall, we find that the influence of the velocity deficit $\delta _d$ is expensive to capture and the systematic study of its influence on the simulation results is left for future work. For the rest of this work, we will only consider the $\delta _d=1$ case in simulations. Nonetheless, we wish to comment on the reason for this difficulty. As pointed out by Otto et al. (Reference Otto, Rossi and Boeck2013), and seen in figure 1(a), the presence of an interfacial velocity deficit induces the presence of a minimum in the velocity profile, in our case on the liquid side, around which the shear changes sign. Depending on the value of $\delta _d$ and other injection parameters, this minimum can move very close to the interface. In our case this minimum is located at a distance of $238\,\mathrm {\mu }$m from the interface, hence a distance two times smaller than the vorticity thickness. Reducing $\delta _d$ or $M$ can make this minimum move even closer to the interface (Otto et al. Reference Otto, Rossi and Boeck2013). The position of this minimum and the associated shear are very expensive to capture in simulations, hence the difficulty to capture the effect of $\delta _d$. Unlike the case of convective modes presented in Bagué et al. (Reference Bagué, Fuster, Popinet, Scardovelli and Zaleski2010), we do not observe the presence of extremely sharp variations of the eigenfunctions associated with the mode of instability, as can be seen on figure 6, where we show $|\varPhi |$, the norm of the stream function computed by the stability analysis for case $B1$ and for two values of $\delta _d$: $\delta _d=1$ (dashed line) and $\delta _d=0.5$ (continuous line), as a function of the vertical coordinate. The stream function is related to the vertical velocity perturbation: $\partial \varPhi /\partial x= \tilde {v}(k,\omega,y)\text {e}^{\text {i}(kx-\omega t)}$. No significant difference is observed between both cases, except a slightly higher value of $|\varPhi |$ in the gas stream for $\delta _d=1$. Thus, we think that the difficulty associated with the resolution of the velocity deficit is solely due to the difficulty of resolving the velocity minimum position and the associated shear, rather than to an increase in the eigenfunctions sharpness. Note that this comment only applies to confinement-induced absolute modes which occur at low wavenumbers (Matas Reference Matas2015).

Figure 6.

Absolute value of the stream function eigenfunction for case B1 with $H_g=H_l=1$ cm: Dotted line, $\delta _d=1$ ($k=140+170i$ m$^{-1}$, $\omega =209+30i$ s$^{-1}$); continuous line, $\delta _d=0.5$ ($k=135+195i$ m$^{-1}$, $\omega =272 +35i$ s$^{-1}$); red, liquid side; blue, gas side. The horizontal axis is shifted so that the $y=0$ location corresponds to the interface location.

This mesh convergence study shows that, focusing on the case where no interfacial velocity deficit is introduced, a resolution $n=8$ allows a converged amplitude and a very low error on the frequency value. We will therefore choose this resolution in § 4.3 where we will compare simulation results obtained for all cases presented in table 2. However, the previous mesh convergence study also shows that a resolution $n=2$ leads to an acceptable error on the frequency, as well as an amplitude that follows the self-similar growth model. We will therefore choose $n=2$ in § 5, devoted to the analysis of the effect of confinement on the stability of the flow.

We have thus established the performance and the limits of our numerical methods, which could be useful for the development and validation of future methods. We now switch to the use of the solver to demonstrate convergence between simulations, linear theory and experiments on confinement-induced absolute modes.

4.3. Convergence between simulation, linear theory and experiments

In figure 7(a), we show the evolution of the most amplified frequency obtained numerically, for $n=8$ ($\circ$, blue) and $n=2$ ($\ast$, orange), as a function of the dynamic pressure ratio and for all cases of table 2. We plot the numerical results along with the most amplified frequency obtained with linear stability analysis ($\triangle$, magenta). For all cases of table 2, the stability analysis predicts that the instability is absolute, driven by a resonance of the shear instability within the confined streams (confinement mechanism). We also add for reference the experimental results of Matas et al. (Reference Matas, Marty and Cartellier2011) ($\square$, black).

Figure 7.

(a) Most amplified frequency of the interface height signal as a function of the dynamic pressure ratio for all cases of table 2 with: experimental results (Matas et al. Reference Matas, Marty and Cartellier2011) ($\square$, black); linear stability analysis results ($\triangle$, magenta); simulation results with $n=2$ ($\ast$, orange); simulation results with $n=8$ ($\circ$, blue). (b) Spectrogram of the interface height signal for case $A1$, with $n=8$. The horizontal magenta line represents the predicted wave frequency.

The frequency values obtained by simulations ($n=8$), stability analysis and experiments are summarized in table 3, along with the error level between the different approaches and the values for each case of the temporal growth rate at the pinching point, $\omega _{i,0}$, also called the absolute growth rate.

Table 3.

Comparison between wave frequency values obtained using simulations, stability analysis and experiments, error between the three approaches, and absolute growth rates for each case.

The comparison between stability analysis, simulations and experiments has to be split in two parts. First, the comparison between stability analysis and simulations, which are directly connected as the inlet velocity profile in simulations is the base flow used in stability analysis, presents a maximal level of error of 11.6 %. This is very satisfactory. The remaining difference might be explained by nonlinearities, spatial variations of the velocity profile or by subtle numerical errors.

Secondly, the comparison of the two previous approaches (simulations and stability analysis) with experiments presents a maximum error of 36.5 % on wave frequency value. Many choices in the modelling can explain such discrepancies, e.g. shape of the velocity profile (which includes velocity deficit value, vorticity thicknesses values, use of ${\rm erf}$ functions), turbulence intensity, nonlinearities, 3-D effects, as well as measurements errors. Given the complexity of the phenomenon, we consider that an average error level of 10 % and a maximum error level of less than 40 % between simulations/stability analysis and experiments is satisfactory.

As shown in table 3, we note that the case with the lowest absolute growth rate corresponds to a case where the agreement between the three approaches is extremely satisfactory. An absolute instability growing exponentially in space and time, the previous observation is not surprising as the case with the lowest absolute growth rate may correspond to the case where nonlinear effects will take the longest before they dominate (as will also be seen in the study of wave amplitude). For the other cases, no correlation seems to emerge between $\omega _{i,0}$ and the error between the three approaches, but we note that these five cases are farther from the neutral boundary curve ($\omega _{i,0}=0$ s$^{-1}$).

To go further in the comparison, one could attempt a global stability study, or several local stability analyses, by taking into account the spatial variations of the base flow, as done for example in Juniper et al. (Reference Juniper, Tammisola and Lundell2011) and Biancofiore et al. (Reference Biancofiore, Gallaire and Pasquetti2011), or directly using nonlinear mean flows, as done for example in Oberleithner, Rukes & Soria (Reference Oberleithner, Rukes and Soria2014). This is something that we have attempted, but we faced difficulties in fitting velocity profiles obtained from simulations by velocity profiles that would still be in the family of velocity profiles described by (2.1)–(2.3). We therefore leave this approach for future work.

This favourable comparison between all approaches is a significant result because it shows that, while the level of error may vary between cases, we do observe convergence between simulation, linear stability analysis and experiment on cases where the instability is triggered by confinement.

The inclusion of the finite thickness of liquid and gas streams in the stability analysis has been shown to allow the convergence between experiments and linear stability analysis (Matas Reference Matas2015; Matas et al. Reference Matas, Delon and Cartellier2018); here we confirm for the first time its impact through simulations. Our work is, therefore, highly complementary to the one presented in Fuster et al. (Reference Fuster, Matas, Marty, Popinet, Hoepffner, Cartellier and Zaleski2013), which showed convergence between simulation, linear stability analysis and experiments on surface tension-induced absolute instabilities. In § 5, we will give further evidence of the destabilizing influence of confinement.

In figure 7(b), we show a spectrogram of the interface height signal for case A1 with $n=8$. One can see that the most amplified frequency is not perfectly constant with downstream distance, but that it matches the theoretical value for more than half of the spatial locations. We can thus question the appearance of such small spatial variations of the most amplified frequency. A possibility would be to question the local or global nature of the mode of instability (Huerre & Monkewitz Reference Huerre and Monkewitz1990). The frequencies predicted by linear stability analysis correspond to the stability property of the inlet velocity profile, but do not give information of whether the predicted absolute mode is local or global. However, it is observed experimentally that the most amplified frequency does not vary with downstream distance (Matas et al. Reference Matas, Marty and Cartellier2011), i.e. that the absolute mode is global. We observed in our simulations that those spatial frequency variations are increasing with resolution. Therefore, we postulate instead that this is an effect of the 2-D nature of the simulation, and particularly of the 2-D nature of the turbulence. Note that a decrease of frequency with downstream distance has already been observed for similar 2-D configurations with moderate density ratios, and was attributed to vortex pairing (Fuster et al. Reference Fuster, Bagué, Boeck, Le Moyne, Leboissetier, Popinet, Ray, Scardovelli and Zaleski2009) or wave merging (Valluri et al. Reference Valluri, Naraigh, Ding and Spelt2010). A way to improve the understanding of spatial frequency variations would be to perform a mesh convergence in three dimensions under experimental conditions, but this task is made difficult by its very high computational cost.

In figure 8, we present a comparison of the wave amplitude evolution for the cases of table 2. In figure 8(a), one can see that for all cases studied, the amplitude follows a similar evolution as the one described in § 4.1. The main differences in terms of amplitude are found for $x/\delta _g \leq 50$, where we see that, for the range of injection conditions under study, the amplitudes are sorted with the dynamic pressure ratio, as detailed by the insert in figure 8(a). Interestingly, one can also note that the absolute growth rate is also increasing with $M$, see table 3. This relation between downstream wave amplitude evolution and absolute growth rate is not surprising as an absolute instability is growing exponentially in space and time. This effect of the dynamic pressure ratio on the near injector amplitude is coherent with the creation of a potential liquid cone that shortens with $M$ (Raynal Reference Raynal1997). Downstream of this zone, the amplitudes are closer to each other, but clearly follow a linear growth with downstream distance, as predicted by (4.2). To further examine this zone of linear growth, we present in figures 8(b) and 8(c) the evolution of the amplitude at fixed gas Reynolds number, $Re_g=14\,667$ and $Re_g=18\,000$, respectively, and for the three values of the liquid Reynolds number used in this study. In both of those figures, we add the self-similar growth model, (4.2), for comparison. The agreement with simulation results is very good. Remarkably, $C_0$ is constant for all plots, and the differences in the slopes are only due to the change in injection conditions, showing a clear self-similarity of the wave amplitude. This was not the case in the experiments of Marty (Reference Marty2015), where $C_0$ had to be adjusted to fit the experimental results. This is consistent with our observation of § 4.2: when the wave speed changes, some corrections have to be applied to the velocity scale used in (4.2). Since all the cases simulated in this section correspond to confinement-induced absolute modes, the wave velocity is close to $U_D$, which is in agreement with the inviscid mechanism driving the instability (Matas Reference Matas2015), therefore, no correction is needed in the self-similar growth model.

Figure 8.

Comparison of wave amplitude for all cases of table 2, for $n=8$ with A1 ($\circ$, blue); A2 ($\circ$, red); A3 ($\circ$, magenta); B1 ($\circ$, orange); B2 ($\circ$, green), B3 ($\circ$, black). (a) evolution of the amplitude for all 6 cases. The insert shows the same data with a zoom in on the inlet area. (b) Evolution of the amplitude for $Re_g=14\,667$. (c) Evolution of the amplitude for $Re_g=18\,000$. In (b,c), the solid lines correspond to the self-similar growth model, (4.2), with the same colour code as the symbols and $C_0=0.21$. (d) Evolution of the position of entry in the nonlinear regime as a function of the dynamic pressure ratio for all cases of table 2 ($\square$, orange). The dashed black line corresponds to a $1/M$ scaling. The insert shows the same data plotted as a function of the gas Reynolds number with: $Re_l=2600$ ($\square$); $Re_l=3700$ ($\times$); $Re_l=5000$ ($\triangle$). Note that in (a–c), only one data point in every six is shown to enhance readability.

In figure 8(d) we present the position of entry in the self-similar region as a function of the dynamic pressure ratio. For the range of injection conditions studied here, one can see that the position of entry in the self-similar zone appears to scale as $1/M$. We also plot this position as a function of the gas Reynolds number, which shows that the impact of the liquid Reynolds number increases when the gas Reynolds number is reduced. Note that, in figure 5(d), one can observe that the decrease of $\delta _d$ seems to induce an earlier entry in the nonlinear growth region.

In summary, we have demonstrated the agreement between linear theory, simulations and experiments on wave frequency value in case of confinement modes, which validates the role of confinement in the stability of an air–water mixing layer. We have also shown the consistency between simulations and experiments on wave growth, as well as a verification of the self-similar model for nonlinear wave growth of Hoepffner et al. (Reference Hoepffner, Blumenthal and Zaleski2011).

5.1. Convective/absolute transition

As already said in the introduction, convectively unstable flows behave as noise amplifiers, i.e. the response of the system depends on the injected perturbations, whereas absolutely unstable flows behave as oscillators, i.e. the response of the system does not depend on the injected perturbations (Huerre & Monkewitz Reference Huerre and Monkewitz1990).

In order to apply those notions to our study, we first focus on case B2, with $H_g=H_l=H$, i.e. a symmetric confinement. We study two values of $H$, namely $H=1$ cm and $H=0.25$ cm. We show in figure 9 that the flow is predicted by linear stability analysis to be absolutely unstable for $H=1$ cm ($\bullet$) and convectively unstable for $H=0.25$ cm ($*$, blue). In the case of the absolute mode, a pinching occurs at low wavenumber (and positive $\omega _i$) between confinement and shear branches. For the convective mode, the most unstable mode of the shear branch, corresponding to the largest $|k_i|$, is located at a higher wavenumber. For $H=0.25$ cm ($*$, blue), the confinement branch can be seen close to the imaginary axis but no pinching occurs with the shear branch.

Figure 9.

Spatial branches for case B2 with: $H=1$ cm ($\bullet$) and $H=0.25$ cm ($*$, blue).

In figure 10, we examine the spectrograms of the instability for these two values of $H$ and different injected perturbations. We choose a uniformly random, low-amplitude perturbation (figure 10a,b), as in the previous section, and a deterministic perturbation $f'=10^{-3} \sin (2{\rm \pi} f_{pert} t)$ (figure 10c–f), with $f_{pert}$, the perturbation frequency to be specified later.

Figure 10.

Noise amplifier versus oscillator behaviour; case B2, $n=2$. Spectral resolution is $0.56$ Hz for all spectrograms. Here (a,c,e) $H=H_g=H_l=1$ cm and (b,d, f) $H=H_g=H_l=0.25$ cm. (a,b) Random noise of amplitude $10^{-5}$. (c–f) Deterministic noise of amplitude $10^{-3}$ and $f_{pert}={15;40;60;80}$ Hz for the cases (c–f), respectively. The colour scale corresponds to the normalized spectral power.

One can see that when the mode is predicted to be absolute ($H=1$ cm, figure 10a,c,e), the spectrograms are barely modified by the injected perturbations. The injected deterministic perturbation can be seen in figures 10(c) and 10(e), only very close to the inlet ($x/\delta _g \leq 5$). Farther downstream, the flow has its own dynamics, and the most amplified frequency and its variations are not affected by the injected perturbation.

On the other hand, when the mode is predicted to be convective ($H=0.25$ cm, figure 10b,d, f), the spectrograms look very different depending on the injected noise. In the case of a random, low-amplitude perturbation, shown in figure 10(b), the most amplified frequency is more difficult to extract as it varies spatially. This behaviour had already been seen in Fuster et al. (Reference Fuster, Bagué, Boeck, Le Moyne, Leboissetier, Popinet, Ray, Scardovelli and Zaleski2009), where it is also said that for the case of convectively unstable flows, the comparison between local linear stability analysis and simulations has to be restricted to the linear growth region, i.e. where the linear stability analysis makes sense. In this area, around $x/\delta _g = 25$, the most amplified frequency is of $73$ Hz. Farther downstream, the most amplified frequency is progressively decreasing. For the cases with a deterministic perturbation injection, shown in figures 10(d) and 10( f), one can see that the forcing frequency is the most amplified in more than half of the domain, showing that the stability of the flow does actually depend on the injected perturbation.

Another way to illustrate that we capture a convective/absolute transition driven by confinement is by constructing a spatiotemporal diagram of the interface height, also called a Hovmöller diagram, as used in Odier et al. (Reference Odier, Balarac, Corre and Moureau2015) and Odier, Balarac & Corre (Reference Odier, Balarac and Corre2018). As shown in figure 9, in the confinement-induced absolute mode, the pinching point occurs at a low wavenumber, whereas in the case of convective mode, the most amplified wavenumber is much higher. This difference can also be seen through the Hovmöller diagram of the instability. In figure 11, we show, for case B2 and a random inlet forcing, two spatiotemporal diagrams of the interface height for the same time window and two values of $H$: $H=1$ cm (figure 11a) and $H=0.25$ cm (figure 11b). One can see that the main difference between both plots is in the wavelength. The waves are much shorter for $H=0.25$ cm than for $H=1$ cm which is consistent with the much larger predicted wavenumber for these conditions. For $H=1$ cm we have $\lambda \approx 4H$, whereas for $H=0.25$ cm we have $\lambda \approx 0.25H$.

Figure 11.

Spatiotemporal diagrams of the interface height. Case B2, $n=2$ with a low amplitude random forcing at the inlet. Here (a) $H=H_g=H_l=1$ cm; (b) $H=H_g=H_l=0.25$ cm. The magenta line corresponds to a wave speed equal to $U_D$, i.e. $U_{\phi }=U_D=1.30\,\textrm {m}\,\textrm {s}^{-1}$.

According to the energy budget presented in Matas (Reference Matas2015), the absolute mode induced by confinement is driven by inviscid stresses, whereas the convective mode is fed by the work of tangential viscous stresses on the interface. This difference can also be seen by examining the effect of $H$ on the wave speed. In figure 11, we add on both plots a magenta line corresponding to a wave speed equal to $U_D$. For $H=1$ cm (figure 11a), this value is in agreement with the numerical result, which may be interpreted as an argument in favour of an inviscid mechanism, as suggested by Matas (Reference Matas2015). On the other hand, for $H=0.25$ cm (figure 11b), the agreement between the numerical results and $U_D$ is less favourable.

These three previous arguments (noise amplifier versus oscillator behaviour; wavenumber value; wave speed) demonstrate that the transition from the convective to the absolute mode is indeed induced by confinement, as predicted by linear stability analysis, and is fully captured numerically. A large wavelength of the order of the injector size and a wave speed close to Dimotakis’ may be good indicators of the presence of confinement modes.

5.2. Effect of confinement on the stability of the flow

In the previous section we confirmed the transition between convective and absolute modes due to flow confinement. We now study the impact of confinement, symmetric or not, on the characteristics of the instability, i.e. frequency and amplitude. We remind the reader that we focus on cases A1 and B2, where we only vary the thicknesses of the gas and liquid streams, $H_g$ and $H_l$, respectively, while keeping all other parameters constant. We first start by considering the effect of $H$ for a symmetric confinement, $H_g=H_l=H$, on the most amplified frequency, shown in figure 12.

Figure 12.

Case of a symmetric confinement: effect of $H$ on the most amplified frequency ($n=2$). (a) Simulation results for case A1 ($\circ$, black), (b) simulation results for case B2 ($\circ$, blue). In both plots the triangles correspond to the frequency predicted by linear stability, with an open symbol for a convective mode and a filled symbol for an absolute mode due to confinement. Spectral resolution is of 0.56 Hz for all simulation results. Dashed lines correspond to a $1/H$ scaling.

First, one can observe that the overall agreement between linear stability analysis and simulations is very satisfactory, in convective as well as in absolute regimes.

One can observe that, when reducing $H$, case A1 ($M=2.32$) presents a transition from absolute to convective regimes for $H$ between $0.75$ cm and $1$ cm. On the other hand, case B2 ($M=6.39$) presents the same behaviour for $H$ between $0.25$ cm and $0.35$ cm. These results suggest that the transition from convective to confinement-induced absolute instabilities does depend on the injection condition, i.e. the value of $M$, as well as on the confining geometry. In other words, when increasing $M$, the flow may be absolutely unstable for a larger range of stream thicknesses.

In figure 13, we show the evolution of the most amplified frequency in the case of non-symmetric streams, i.e. we keep one stream thickness constant while varying the other one. In figure 8 of Matas (Reference Matas2015), it is shown through stability analysis that both thicknesses have symmetric effects on the location of the pinching point, which may be surprising as a flow with a thin gas stream and a large liquid stream could seem very different from a flow with a thick gas stream and a thin liquid stream. This can be related to the symmetry found in the dispersion relation for the sinuous mode of a double shear layer in Juniper (Reference Juniper2006). Here we do confirm this symmetric behaviour through simulations: one can see in figure 13 that in most cases, the results for $H_g=1$ cm, $H_l=H$ (symbols $\circ$ and $\circ$, blue) and $H_g=H$, $H_l=1$ cm (symbols $\square$ and $\square$, blue) are in agreement. For $H \leq 1$ cm, the most amplified frequencies are very similar to the ones presented in figure 12 for the case of symmetric streams, i.e. $H_g=H_l$. Above $H=1$ cm, and in case of absolute mode, the most amplified frequency appears to be constant with $H$, which is not the case for a symmetric confinement (see figure 12).

Figure 13.

Case of a non-symmetric confinement: effect of $H$ on the most amplified frequency ($n=2$). (a) Case A1, simulation results with $H_g=1$ cm and $H_l=H$ ($\circ$, black); simulation results with $H_g=H$ and $H_l=1$ cm ($\square$, black). (b) Case B2, simulation results with $H_g=1$ cm and $H_l=H$ ($\circ$, blue), simulation results with $H_g=H$ and $H_l=1$ cm ($\square$, blue). In both plots the triangles correspond to the frequency predicted by linear stability for $H_g=1$ cm and $H_l=H$, with an open symbol for a convective mode and a filled symbol for an absolute mode due to confinement. The spectral resolution is of 0.56 Hz for all simulation results. Dashed lines correspond to a $1/\min (H_g,H_l)$ scaling.

For case A1, $H_g=1$ cm and $H_l=2$ cm, as well as for $H_l=1$ cm and $H_g=2$ cm, the mode is predicted to be convective, hence the high value (71 Hz) for the frequency predicted by linear stability analysis, whereas for both non-symmetric cases we find numerically a resonance at low frequency (26 Hz). Compared with the case of a symmetric confinement, the fact that a non-symmetric confinement may favour a transition to a convective instability is explained by a decrease in the absolute growth rate at the pinching point. Physically speaking, the absolute instability driven by confinement being similar to a mechanism of resonance, as described in the introduction, we postulate that this resonance will be stronger when both gas and liquid streams have the same thicknesses, since velocity and pressure perturbations in both streams can contribute to the feedback needed for the resonance to occur. Healey (Reference Healey2009) has also seen that a symmetric confinement could favour absolute instabilities compared with an asymmetric confinement. For this case, the pinching between shear and confinement branches occurs at a frequency of 26 Hz, i.e. the value found numerically, but at an absolute growth rate of $-20$ s$^{-1}$. The predicted instability can therefore not be considered as absolute and the value retained for the prediction is therefore that corresponding to the largest $|k_i|$, namely 71 Hz. Note that we refined the mesh, up to $n=8$, but still captured the absolute mode. This discrepancy between numerics and linear stability analysis in the transition from convective to absolute modes may be due to nonlinearities, spatial variations of the base flow, or even to subtle numerical errors (Cossu & Loiseleux Reference Cossu and Loiseleux1998).

Note again that the linear stability results presented in figure 13 correspond to the case $H_g=1$ cm and $H_l=H$, as the symmetric case gives the same result, which has been carefully checked, see also figure 8 of Matas (Reference Matas2015).

One important observation is that the gas stream thickness, $H_g$, has its own effect on the instability via the triggering of this absolute instability. It does not only affect the instability through its effect on the gas vorticity thickness $\delta _g$, unlike what has been observed in Matas et al. (Reference Matas, Marty and Cartellier2011) based on the results of Ben Rayana (Reference Ben Rayana2007). However, both effects are related in the experiments and controlling one independently of the other may not be technically feasible.

In figure 14, we present the effect of confinement on wave amplitude. We start by studying in figure 14(a) the effect of a reduction of $H$ for a symmetric confinement, from $H=2$ cm ($+$, red) to $H=0.15$ cm ($\square$, magenta). One can see that the reduction of $H$ induces a global decrease of the amplitude in the downstream region. The growth close to the inlet, i.e. $x/\delta _g \le 25$, does not seem to be affected until the amplitude is of the order of the gas vorticity thickness. Farther downstream, the lower the value of $H$, the earlier the amplitude seems to depart from the self-similar growth model, and stabilize up to a constant value increasing with $H$.

Figure 14.

Effect of confinement on the amplitude for case B2 with $n=2$. (a) Effect of the reduction of $H$ for a symmetric confinement ($H_g=H_l=H$): $H=2$ cm (+, red); $H=1$ cm ($\circ$, blue); $H=0.5$ cm ($\circ$, green); $H=0.25$ cm ($\triangle$, black); $H=0.15$ cm ($\square$, magenta). (b) Evidence of the non-symmetric effect of the gas and liquid thicknesses on the amplitude: $H_l=H_g=1$ cm ($\circ$, blue); $H_g=0.15$ cm and $H_l=1$ cm ($\circ$, green); $H_l=0.15$ cm and $H_g=1$ cm ($\triangle$, black); $H_l=H_g=0.15$ cm ($\square$, magenta). In both plots, only one data point in every four is shown to enhance readability. Here (4.1), with $A_0/\delta _g=0.0003$ and $k_i=614.6$ m$^{-1}$ (orange dashed line); (4.2), with $C_0=0.21$ (Orange solid line); (4.2), with (a): $C_0=0.045$, and (b): $C_0=0.032$ (Cyan solid line).

To clarify the effect of each stream thickness on the wave amplitude, we present in figure 14(b) the effect of a reduction of $H_g$ only ($\circ$, green), $H_l$ only ($\triangle$, black) and both ($\square$, magenta). One can see that, although the effect of both stream thicknesses on the frequency is symmetric, the effect on the amplitude is non-symmetric. It seems that the decrease in amplitude observed in figure 14(a) between $H=2$ cm (symbol $+$, red) and $H=0.15$ cm (symbol $\square$, magenta) is mainly due to the decrease of $H_l$ since the same decrease is observed when only $H_l$ is reduced (symbol $\triangle$, black in figure 14b). It also seems that when $H_l$ only is reduced, the growth is slower in the region close to the injector, which agrees with the experimental observations of Marty (Reference Marty2015).

Similarly, the departure from the self-similar model seems to be due to the reduction of $H_g$. Following the derivation of the self-similar model in Hoepffner et al. (Reference Hoepffner, Blumenthal and Zaleski2011), the amplitude increase in the nonlinear area is due to the aerodynamic suction created by the acceleration in the gas stream when deviated by the wave. When the wave amplitude becomes large compared with the gas stream, there is no longer a sufficient acceleration to induce wave growth.

5.3. Scaling law for confinement modes

The destabilizing effect of confinement has been demonstrated and explored by means of numerical simulations and linear stability analysis in previous subsections. In Matas et al. (Reference Matas, Delon and Cartellier2018), a scaling law derived from the linearized equations was proposed for the wave frequency in confinement-induced absolute modes,

(5.1)\begin{equation} f \sim \frac{\sqrt{\dfrac{\rho_g}{\rho_l}\dfrac{\delta_l}{\delta_g}}U_g+U_l}{L}, \end{equation}

where $L$ is a characteristic length related to the confinement. This scaling has recently been found to account very well for the effect of the nozzle size in experiments on a coaxial configuration (Singh et al. Reference Singh, Kourmatzis, Gutteridge and Masri2020). Equation (5.1) predicts a linear evolution of the frequency with the liquid velocity, which has been observed experimentally (Matas et al. Reference Matas, Marty and Cartellier2011; Singh et al. Reference Singh, Kourmatzis, Gutteridge and Masri2020). In the case of large gas velocity, $\delta _g \ll \delta _l$, and for $\delta _g$ varying as $U_g^{-1/2}$, the frequency is proportional to $U_g^{5/4}$, which is close to the scaling predicted by inviscid theory (Raynal Reference Raynal1997; Marmottant & Villermaux Reference Marmottant and Villermaux2004), see § 1. Note that the scaling of (5.1) is only valid for absolute modes due to confinement.

In order to compare our results to the scaling of (5.1), we add lines enhancing the scaling of the frequency with $H$ in absolute regimes in figure 12. In the confinement-induced absolute regime, we find that the most amplified frequency follows a ${\sim }1/H$ scaling in the case of symmetric confinement. For $H_g \ne H_l$, as seen in figure 13, it seems that the frequency scales with the most restrictive confining length, i.e.

(5.2)\begin{equation} f \sim 1/\min(H_g,H_l). \end{equation}

Thus, our results match very well with the scaling of (5.1) and we evidence that $L$ is of the order of the smallest injector size. This is different from what has been observed experimentally by Delon, Cartellier & Matas (Reference Delon, Cartellier and Matas2018) for a coaxial geometry, where the frequency is found to always scale with $1/H_l$. The shear branch being barely affected by the value of $H$, frequency is expected to be independent of the confining size in the convective regime.