2.1 Flow configuration and experimental set-up

Figure 1(a) sketches the configuration under study, where a liquid of density $\unicode[STIX]{x1D70C}$ , kinematic viscosity $\unicode[STIX]{x1D708}$ and surface tension coefficient $\unicode[STIX]{x1D70E}$ is injected at a constant flow rate $Q$ through a circular tube of radius $R$ whose length is sufficiently large to guarantee a fully developed velocity profile at its outlet. The resulting free liquid jet is confined in the axial direction through its impact with the free surface of a reservoir of the same liquid placed at a distance $L$ from the injector outlet. To ensure that the results are independent of the downstream boundary condition, we have also performed several experiments where the jet impacts a solid horizontal plate, obtaining nearly identical results. The surrounding gaseous atmosphere, at pressure $p_{a}$ , has a negligible dynamic effect on the jet due to the smallness of the typical liquid velocities.

Figure 1. (a) Sketch of the flow configuration. (b) Photograph of a steady jet for $\unicode[STIX]{x1D708}=500$ cSt $(\unicode[STIX]{x1D6E4}=8.33)$ , $R=1.75~\text{mm}\;(Bo=1.38)$ , $Q=3.5~\text{ml}~\text{min}^{-1}\;(We=3.05\times 10^{-3})$ and $L=16.6~\text{mm}$ $(L/R=9.47)$ . The white line shows the steady solution of (2.1)–(2.3).

The experiments were performed in the same set-up used by Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013), with the only addition of a vertical positioning stage located below the injection tube and mounting a platform onto which a reservoir filled with the working liquid was placed. The free surface of the reservoir, located at a distance $L$ from the injector outlet, was impinged by the liquid jet as shown in figure 1. The liquid was supplied with a Harvard Apparatus PhD Ultra syringe pump, and the free jet was recorded using a Red Lake Motion Pro X high speed camera. To minimise the effect of ambient noise, both the injector and the camera were installed inside a transparent isolation chamber, and the entire system except the syringe pump was placed on an optical table with a passive vibration damping system. Two stainless steel capillary tubes acquired from Tubca were used as injectors, of inner radii $R=1.5~\text{mm}$ and $R=1.75~\text{mm}$ , that were carefully machined at their tip to ensure that the contact line remained pinned at their inner radii. Two different Newtonian liquids were used in the experiments, both of them PDMS (Polydimethylsiloxane) silicon oils from Sigma-Aldrich, whose properties at $25\,^{\circ }\text{C}$ are $\unicode[STIX]{x1D70C}=970~\text{kg}~\text{m}^{-3}$ , $\unicode[STIX]{x1D70E}=21.1~\text{mN}~\text{m}^{-1}$ and kinematic viscosities $\unicode[STIX]{x1D708}=500$ and $\unicode[STIX]{x1D708}=1000~\text{mm}^{2}~\text{s}^{-1}$ , with corresponding values of the Kapitza number, $\unicode[STIX]{x1D6E4}=3\unicode[STIX]{x1D708}[(\unicode[STIX]{x1D70C}^{3}g)/\unicode[STIX]{x1D70E}^{3}]^{1/4}$ , of $\unicode[STIX]{x1D6E4}=8.33$ and $\unicode[STIX]{x1D6E4}=16.67$ , respectively. In each experimental run, first an injection tube and a silicone oil were selected, and the vertical stage was positioned to provide the desired value of $L$ . The syringe pump was programmed to inject an initial liquid flow rate large enough to provide a jetting regime, and to smoothly decrease until the desired target value of $Q$ was achieved.

2.2 Basic experimental evidence

From the results of the experiments it is deduced that, leaving the presence of coiling apart, there are three possible regimes whose occurrence depends on the values of the control parameters. Let us begin by classifying these regimes following the order in which they are observed when $Q$ decreases smoothly and the rest of the control parameters are fixed.

1. (i) Steady jetting: This regime is illustrated in movie 1 of the supplementary material and in figure 1(b), which shows a photograph of a steady jet of silicone oil with $\unicode[STIX]{x1D708}=500$ cSt injected through a tube of radius $R=1.75$ mm at a constant flow rate $Q=3.5~\text{ml}~\text{min}^{-1}$ that impinges on the free surface of an oil reservoir placed at a distance $L=16.6~\text{mm}$ of the injector. The steady jetting regime is stable if $Q$ is larger than a certain critical value, $Q_{c}$ .

When $Q<Q_{c}$ the jet is unstable due to an oscillatory axisymmetric global mode, leading to self-sustained oscillations of shape and velocity whose amplitude increases with time (Sauter & Buggisch Reference Sauter and Buggisch2005; Rubio-Rubio et al. Reference Rubio-Rubio, Sevilla and Gordillo2013). This linear global mode is also referred to as the breathing mode throughout the paper. Our experiments reveal that two different nonlinear states may take place if $Q<Q_{c}$ :

1. (ii) Oscillatory jetting: When $Q_{b}<Q<Q_{c}$ a limit-cycle state without breakup is achieved, whereby the oscillation amplitude saturates to a certain function of the downstream position (see figure 5). This regime can be observed in movies 2 and 5 of the supplementary material, as well as in figures 3–7. Note that the oscillatory jetting regime is a nonlinear state of the liquid thread, and it should not be confused with the linear breathing mode, which is globally unstable under the conditions where both the oscillatory jetting and dripping regimes take place.

1. (iii) Dripping: When $Q<Q_{b}$ , the amplitude of the oscillations grows until the breakup of the jet takes place, finally leading to a jetting–dripping transition, as illustrated in movies 3 and 4 of the supplementary material, as well as in figure 8.

This evidence calls out for an experimental and numerical characterisation of the functions $Q_{c}(\unicode[STIX]{x1D708},R,L)$ and $Q_{b}(\unicode[STIX]{x1D708},R,L)$ , that were thus obtained in a wide region of parameter space, as reported in § 3.

2.3 One-dimensional model

To model the liquid jet we use the dimensionless leading-order one-dimensional mass and momentum equations (Eggers & Dupont Reference Eggers and Dupont1994; García & Castellanos Reference García and Castellanos1994),

where the dependent variables $u$ and $r_{j}$ are the liquid velocity and jet radius respectively, $z$ is the axial coordinate, $t$ is the time and ${\mathcal{C}}$ is the interfacial curvature. It is worth pointing out that the full expression for the curvature needs to be retained for the model to provide good quantitative predictions of the stability in the case of large injector diameters, for which the jet experiences a strong gravitational stretching close to the neutral conditions and, correspondingly, the stabilising effect of the axial curvature must be taken into account (Rubio-Rubio et al. Reference Rubio-Rubio, Sevilla and Gordillo2013). The variables in (2.1)–(2.3) have been made dimensionless using the liquid density $\unicode[STIX]{x1D70C}$ , the capillary length $l_{\unicode[STIX]{x1D70E}}=(\unicode[STIX]{x1D70E}/\unicode[STIX]{x1D70C}g)^{1/2}$ and the associated characteristic velocity $\sqrt{gl_{\unicode[STIX]{x1D70E}}}$  (Senchenko & Bohr Reference Senchenko and Bohr2005). Thus, the system (2.1)–(2.3) only depends on $\unicode[STIX]{x1D6E4}$ , which is constant for a given liquid and a fixed value of $g$ . The solution depends on the constant flow rate, $Q$ , and on the injector radius, $R$ , through the dimensionless versions of the boundary conditions at $z=0:r_{j}=R/l_{\unicode[STIX]{x1D70E}}=Bo^{1/2}$ and $u=U/\sqrt{gl_{\unicode[STIX]{x1D70E}}}=We^{1/2}Bo^{-1/4}$ , where $Bo=\unicode[STIX]{x1D70C}gR^{2}/\unicode[STIX]{x1D70E}$ is the Bond number, $We=\unicode[STIX]{x1D70C}U^{2}R/\unicode[STIX]{x1D70E}$ is the Weber number and $U=Q/(\unicode[STIX]{x03C0}R^{2})$ the mean velocity at the nozzle exit. In addition, a Dirichlet boundary condition is imposed for the velocity at the downstream end of the domain, $z=L/l_{\unicode[STIX]{x1D70E}}$ : $u=u_{out}$ , as a crude but simple way to represent the impingement of the jet on the reservoir. Specifically, the value of $u_{out}$ is small enough to properly describe the impact region, since the liquid bath is at rest far from the jet. This simple method provides fairly good results, as illustrated in figure 1(b) for the particular case of a steady jet. In addition, we have carefully checked that the results barely depend on the value chosen for the outflow velocity, provided that $u_{out}\lesssim 0.1$ . In the case of relatively long jets with $L/R\gtrsim 20$ , the value of $u_{out}$ does not affect the results at all. Indeed, in this limit the impact region is very small compared to the length of the jet, and the outflow boundary condition affects neither the base flow nor its linear and nonlinear stability. In particular, these cases can be easily modelled either by imposing a Neumann outlet boundary condition, or even without imposing any outlet boundary condition at all (Rubio-Rubio et al. Reference Rubio-Rubio, Sevilla and Gordillo2013). The mathematical model is governed by four control parameters, namely $\unicode[STIX]{x1D6E4}$ , $Bo$ and $We$ and the dimensionless length of the jet, $L/R=L/l_{\unicode[STIX]{x1D70E}}Bo^{-1/2}$ , hereafter scaled with the injector radius for convenience.

3.1 Linear stability analysis: the role of axial confinement

Let us begin by extending the linear stability results of Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013) to account for the effect of $L/R$ on the critical Weber number, $We_{c}$ , below which the jet becomes globally unstable. Since the methodology is identical to that presented in Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013), except for the treatment of the downstream boundary condition, the details of the formulation are provided in appendix A. The procedure consists of linearising equations (2.1)–(2.3) around a given steady basic state according to the following decomposition in temporal normal modes,

where $\unicode[STIX]{x1D716}\ll 1$ accounts for the smallness of the perturbation around the base flow $[r_{j0}(z),u_{0}(z)]$ , $\unicode[STIX]{x1D714}$ is the complex eigenfrequency and $r_{j1}(z)$ and $u_{1}(z)$ are the corresponding eigenfunctions. The values of $\unicode[STIX]{x1D714}_{r}=\text{Re}(\unicode[STIX]{x1D714})$ and $\unicode[STIX]{x1D714}_{i}=\text{Im}(\unicode[STIX]{x1D714})$ represent the growth rate and the angular frequency of a normal mode, respectively. Accordingly, the jet will be stable if $\max (\unicode[STIX]{x1D714}_{r})<0$ , and unstable if $\max (\unicode[STIX]{x1D714}_{r})>0$ . In the latter case, it is important to emphasise that the linearised description given by (3.1)–(3.2) is only valid during the initial stages of growth of small disturbances, so that the nonlinear terms are negligible in (2.1)–(2.3). It is also worth pointing out that the decomposition (3.1)–(3.2) does not rely on a local stability analysis, which would involve the usual expansion in normal modes of wavepacket form $\text{e}^{\unicode[STIX]{x1D714}t+\text{i}kz}$ . Thus, instead of searching for a dispersion relation $D(\unicode[STIX]{x1D714},k)=0$ , here no a priori assumption is made about the shape of the eigenfunctions, whose spatial structure is obtained as part of the solution.

The basic flow, which satisfies the steady version of (2.1)–(2.3), was calculated using the procedure detailed in § A.1. Once the base flow has been found, the critical conditions for the transition from a globally stable to a globally unstable jet are easily determined as a function of the governing parameters, $(\unicode[STIX]{x1D6E4},Bo,We,L/R)$ , by solving the linear stability problem as explained in § A.2. If the flow is stable, $\max (\unicode[STIX]{x1D714}_{r})<0$ , small disturbances are damped and thus a steady flow is expected, like that shown in figure 1(b) and in movie 1 of the supplementary material. In contrast, if $\max (\unicode[STIX]{x1D714}_{r})>0$ the jet is globally unstable and the development of spontaneous oscillations of increasing amplitude is predicted. In the latter case, two different scenarios emerge according to the experimental evidence described in § 2.2: either the oscillations saturate to a limit cycle without breakup (see figures 3–7, and movies 2 and 5 of the supplementary material), or their growth leads to the pinch-off of the liquid column, eventually leading to a dripping regime (see figure 8 and movies 3 and 4 of the supplementary material). Since the oscillation amplitude increases in the downstream direction (see figure 5) it can be anticipated that the confinement parameter, $L/R$ , will have a strong influence not only on the neutral stability curves, but also on the selection of the final nonlinear regime (see § 3.2).

Figure 2. (a) Critical Weber number, $We_{c}$ , and (b) corresponding angular frequency, $\unicode[STIX]{x1D714}_{c,i}$ , as functions of $L/R$ , for $(\unicode[STIX]{x1D6E4},Bo)=(16.67,1.00)$ (solid lines and filled symbols) and $(\unicode[STIX]{x1D6E4},Bo)=(8.33,1.38)$ (dashed lines and open symbols). The lines are the neutral curves obtained with the linear stability analysis. The symbols represent the experimental results for the maximum (▿) and minimum (▵) values of $We_{c}$ , and for the critical frequency (○). The vertical thin lines in panel (a) are the critical lengths given by the hydrostatic limit, $We_{c}\rightarrow 0$ (Kovitz Reference Kovitz1975; Benilov & Cummins Reference Benilov and Cummins2013).

Both the experiments and the stability analysis reveal that the eigenvalue spectrum is strongly affected by the axial confinement for sufficiently small values of $L/R$ . Consequently, $We_{c}$ and $\unicode[STIX]{x1D714}_{c,i}$ are certain functions of $L/R$ which can be easily computed with the methodology described in §§ A.2 and B.1. The main result is summarised in figure 2, which shows the good agreement of the experiments (symbols) with the prediction of the linear stability analysis (lines), especially for values of $L/R\lesssim 10$ . Figure 2 reveals that both $We_{c}$ and $\unicode[STIX]{x1D714}_{c,i}$ decrease as $L/R$ decreases, indicating that axial confinement stabilises the jet and reduces the critical frequency of the self-sustained oscillations. The symbols ▵, ▿ in figure 2(a) define the experimentally determined range of neutral stability, due to the fact that the flow rate decreased in smooth ramps programmed with the syringe pump as described in § 2.1. For sufficiently long jets, $L/R\gg 1$ , the stability properties become independent of $L/R$ , reaching the limit already studied by Sauter & Buggisch (Reference Sauter and Buggisch2005) and Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013). In the opposite limit of strongly confined jets, figure 2(a) suggests the existence of a vertical asymptote, $We_{c}\rightarrow 0$ , and a corresponding critical value of the jet length, $L_{c}/R$ . Indeed, the vertical lines plotted in figure 2(a) are the maximum lengths for which a static axisymmetric liquid bridge between a solid rod and an infinite pool is stable (Kovitz Reference Kovitz1975; Benilov & Cummins Reference Benilov and Cummins2013), namely $(Bo,L_{c}/R)\simeq (1,1.31)$ and $(1.38,1.19)$ . Note from figure 2 that these values of $L_{c}/R$ are consistent with our results in the hydrostatic limit, $We_{c}\rightarrow 0$ . It is therefore concluded that axial confinement stabilises the liquid thread, and that the jet length below which such effect is noticeable depends on $\unicode[STIX]{x1D6E4}$ and $Bo$ , having values $L/R\lesssim 15$ for $(\unicode[STIX]{x1D6E4},Bo)=(8.33,1.38)$ and $L/R\lesssim 23$ for $(\unicode[STIX]{x1D6E4},Bo)=(16.67,1.00)$ , as can be deduced from figure 2. In contrast, the hydrostatic limit reached as $L/R$ decreases is independent of the parameter $\unicode[STIX]{x1D6E4}$ , which incorporates the liquid viscosity, and is only a function of $Bo$ . Figure 2 also reveals that the quantitative agreement between experiments and theory deteriorates for $L/R\gtrsim 10$ , probably due to the limitations of the one-dimensional model. In particular, as emphasised by Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013), the model does not contemplate the viscous relaxation from the parabolic velocity profile at the injector outlet to the uniform velocity profile achieved downstream. In contrast, the good quantitative agreement found for $L/R\lesssim 10$ may well be due to the fact that the hydrostatic limit, $We_{c}\rightarrow 0$ , is described exactly by the theoretical model thanks to the use of the full expression for the interfacial curvature in (2.3).

3.2 Nonlinear stability

The present section is devoted to address the influence of the control parameters on the selection of the final jet state under globally unstable conditions, $Q<Q_{c}$ . Although the linear stability analysis presented in § 3.1 provides values of $Q_{c}$ in fairly good agreement with experiments, it cannot predict the final state of the jet at large times, which is determined by nonlinear effects. Therefore, in addition to the experiments, we have conducted numerical simulations of (2.1)–(2.3), which were integrated by means of the simple and efficient method explained in §§ B.1 and B.2. In particular, the latter approach allows us to numerically compute the nonlinear saturation process that takes place in the oscillatory jetting regime, and to determine the breakup flow rate, $Q_{b}$ , as a function of the governing parameters $\unicode[STIX]{x1D708}$ , $R$ and $L$ .

Let us first describe the main characteristics of the oscillatory jetting regime, which takes place for $Q_{b}<Q<Q_{c}$ or, in dimensionless terms, for $We_{b}<We<We_{c}$ . For clarity, hereinafter the values of the variables provided by experiments, linear stability analysis and numerical simulations will be denoted using the superscripts $(\text{})^{exp}$ , $(\text{})^{lsa}$ and $(\text{})^{num}$ , respectively.

Figure 3. Logarithm of the radius oscillation amplitude at $z=z^{\ast }=0.73L$ as a function of time extracted from a numerical simulation of (2.1)–(2.3) starting from a very small disturbance superimposed on a slightly unstable steady solution, namely $\unicode[STIX]{x1D6E4}=16.67$ , $Bo=1$ , $L/R=23.26$ and $We^{num}=6.8\times 10^{-3}<We_{c}^{num}$ . After an initial stage of exponential growth during the linear regime, $t\lesssim 800$ , the amplitude saturates to a constant value for $t\gtrsim 2000$ due to nonlinear effects. The inset shows the power spectral density of the saturated oscillations, where the most energetic frequency and its two leading harmonics can be appreciated.

Figure 3 shows a numerical example of the oscillation amplitude growth induced by a very small initial disturbance superimposed on the steady state solution under globally unstable conditions, namely $\unicode[STIX]{x1D6E4}=16.67$ , $Bo=1$ , $L/R=23.26$ and $We^{num}=6.80\times 10^{-3}<We_{c}^{num}=We_{c}^{lsa}=7.88\times 10^{-3}$ . In figure 3 the logarithm of the maxima of a pointwise measure of the radius oscillation, $A_{r_{j}}=r_{j}(z^{\ast },t)-r_{j0}(z^{\ast })$ , where $z^{\ast }=0.73L$ , is plotted as a function of time. The choice of the value $z^{\ast }=0.73L$ to measure the jet oscillations was motivated by the fact that their amplitude is sufficiently large at this distance from the injector to allow a precise measurement. In addition, in cases where coiling takes place, like that shown in figure 5(b), the coiling amplitude at this point is small enough for its influence on the results to be negligible. Three different stages can be clearly identified in figure 3. First, an initial linear growth regime for $t\lesssim 800$ , where the small disturbance increases exponentially with time with a growth rate given by $\unicode[STIX]{x1D714}_{r}^{num}=3.28\times 10^{-3}$ and $\unicode[STIX]{x1D714}_{r}^{lsa}=3.23\times 10^{-3}$ according to the numerical simulation and the linear stability analysis presented in § 3.1, respectively. During the second stage, $800\lesssim t\lesssim 2000$ , a transition regime takes place where nonlinear effects become important and moderate the exponential growth. Third, for $t\gtrsim 2000$ , the oscillation amplitude saturates to a certain constant. The inset of figure 3 represents the power spectral density (PSD) of the saturated signal, clearly showing the dominant angular frequency and its two leading harmonics. The value of the dominant frequency extracted from the numerical PSD is $\unicode[STIX]{x1D714}_{i}^{num}=0.259$ , to be compared with the corresponding value of the leading eigenmode computed with the linear stability analysis of § 3.1, $\unicode[STIX]{x1D714}_{i}^{lsa}=0.259$ .

Figure 4. (a) Photographs of the liquid jet during a period of self-sustained oscillations for $\unicode[STIX]{x1D708}=500$ cSt ( $\unicode[STIX]{x1D6E4}=8.33$ ),  $R=1.75~\text{mm}\;(Bo=1.38)$ , $L=28~\text{mm}\;(L/R=16)$ and $Q^{exp}=3.5~\text{ml}~\text{min}^{-1}\;(We^{exp}=2.96\times 10^{-3})$ , smaller than the experimental critical flow rate, $Q_{c}^{exp}=3.8~\text{ml}~\text{min}^{-1}\;(We_{c}^{exp}=3.49\times 10^{-3})$ . The time interval between photographs is 38 ms, and the oscillation frequency is 3.33 Hz ( $\unicode[STIX]{x1D714}_{i}^{exp}=0.258$ ). (b) Numerically computed period of self-sustained oscillations for $\unicode[STIX]{x1D6E4}=8.33$ , $Bo=1.38$ , $L/R=16$ and $We^{num}=4.43\times 10^{-3}$ , smaller than the critical flow rate provided by the one-dimensional model, $We_{c}^{lsa}=We_{c}^{num}=5.03\times 10^{-3}$ . The time interval is $37$ ms and the oscillation frequency is 3.83 Hz ( $\unicode[STIX]{x1D714}_{i}^{num}=0.297$ ). Note that similar distances to the critical point have been chosen in the experiment and in the numerical simulation, namely $We_{c}^{exp}-We^{exp}\simeq 5.3\times 10^{-4}$ and $We_{c}^{num}-We^{num}\simeq 6\times 10^{-4}$ in (a,b), respectively.

The typical limit-cycle behaviour observed at large times in the oscillatory jetting regime is illustrated in figure 4. Specifically, figure 4(a) shows a sequence of photographs captured during one period of self-sustained oscillations, and figure 4(b) displays the corresponding numerically computed interface profiles under similar conditions. The values of all the dimensionless parameters except $We$ are the same in figure 4(a,b), namely $\unicode[STIX]{x1D6E4}=8.33$ , $Bo=1.38$ , $L/R=16$ , for which the critical Weber numbers for the onset of global instability are $We_{c}^{exp}=3.49\times 10^{-3}$ and $We_{c}^{lsa}=We_{c}^{num}=5.03\times 10^{-3}$ according to the experiments and to the prediction given by the linear stability analysis and the numerical integration of (2.1)–(2.3), respectively. Due to the difference in the experimental and numerical values of $We_{c}$ , which can also be observed in figure 2(a), we decided to choose values of $We^{exp}$ and $We^{num}$ in figure 4 such that the distance to the corresponding critical Weber numbers was similar, namely $We_{c}^{exp}-We^{exp}\simeq 5.3\times 10^{-4}$ in the experiment of figure 4(a) and $We_{c}^{num}-We^{num}\simeq 6\times 10^{-4}$ in the numerical simulation of figure 4(b). In view of the results shown in figure 4 it is clear that the simple one-dimensional model used in the present work is able to qualitatively reproduce the main features of the oscillatory jetting regime including not only its frequency, but also the spatial structure of the limit cycle.

Figure 5. Upper and lower envelopes of the oscillations around the basic steady state, $r_{j}-r_{j0}$ , after the saturation to a limit cycle, obtained from experiments (dotted lines) and numerical simulations of (2.1)–(2.3) (solid lines). Numerical jet shapes at several instants during one oscillation period are also plotted as thin solid lines bounded by their corresponding envelope. The inset plots show the difference between the maxima and minima of $r_{j}$ as a function of $z$ . (a) Configuration without coiling, illustrated by the inset photograph: $\unicode[STIX]{x1D6E4}=8.33$ , $Bo=1.38$ , $L/R=9.36$ and values of $We<We_{c}$ close to the experimental and numerical neutral conditions, namely $We^{exp}=2.96\times 10^{-3}$ , for which $\unicode[STIX]{x1D714}_{i}^{exp}=0.269$ , and $We^{num}=3.75\times 10^{-3}$ , for which $\unicode[STIX]{x1D714}_{i}^{num}=0.296$ . (b) Configuration with coiling, illustrated by the inset photograph: $\unicode[STIX]{x1D6E4}=16.67$ , $Bo=1$ , $L/R=23.26$ and values of $We^{exp}=3.93\times 10^{-3}$ , for which $\unicode[STIX]{x1D714}_{i}^{exp}=0.215$ , and $We^{num}=6.42\times 10^{-3}$ , for which $\unicode[STIX]{x1D714}_{i}^{num}=0.256$ . The vertical dashed line marks the upper limit imposed on the value of $z$ in the post-processing of the experiment.

A quantitative comparison between the typical experimental and numerical behaviour in the oscillatory jetting regime is provided in figures 5–7. In particular, figure 5 shows the numerically computed saturated oscillations of the jet radius around the basic steady state, $r_{j}(z,t)-r_{j0}(z)$ as $t\rightarrow \infty$ (thin solid lines). The thick solid and dashed lines are, respectively, the upper and lower envelopes of the numerical and experimental oscillations, obtained as the maximum and minimum values of $r_{j}-r_{j0}$ over time at each value of $z$ . The insets display the amplitude of the oscillations as the difference between the upper and the lower envelopes, $\max (r_{j})-\min (r_{j})$ . The oscillation amplitude is seen to increase monotonically downstream until the impact region is reached, explaining why the liquid jets studied herein always break up close to the free surface of the reservoir, as observed in figure 8 and in movies 3 and 4 of the supplementary material. A particularly interesting aspect of the flow is illustrated in figure 5(b), which shows that the self-sustained oscillations may coexist with the phenomenon of coiling (Ribe et al. Reference Ribe, Habibi and Bonn2012). Indeed, our experiments have revealed new regimes of steady and oscillatory coiling as a function of the parameters of the problem, respectively associated with the steady jetting and with the oscillatory jetting regimes. In the latter case, the coiling is unsteady and its frequency varies enslaved to that of the axisymmetric breathing mode. Another scenario that we have observed is the disappearance of coiling due to the increase of the oscillation amplitude as the value of $We$ decreases or the value of $L/R$ increases (see figure 9 a,c). Although these features cannot be predicted using an axisymmetric model, their effect on the saturation amplitude upstream of the impact region is relatively small, as deduced from the results of figure 5(b). In fact, it is deduced from figure 5 that the quantitative agreement between the experiments and the numerical integration of the one-dimensional model equations (2.1)–(2.3) is fairly good, provided that the values of $We$ are chosen such that $We_{c}-We$ have similar values in the experiments and numerical simulations.

Figure 6. Axial coordinate of the point of maximum radius of the liquid bulge, $z_{max}$ , as a function of time during one cycle of the self-sustained oscillatory jetting state of figure 5(b), obtained from experiments (symbols) and numerical simulations (solid line). The dashed line shows the free-fall law $t^{2}/2$ .

Another interesting aspect of the oscillatory jetting regime is the periodic formation of a liquid bulge which falls downstream during each oscillation period, as deduced from figures 4 and 5, and from movies 2 and 5 of the supplementary material. Figure 6 shows the axial position of the point of maximum radius inside the liquid bulge, $z_{max}$ , as a function of time, under the conditions of figure 5(b). The symbols and the solid line represent the experimental and numerical results, respectively, while the dashed line is a plot of the free-fall law which, in dimensionless terms, is the function $t^{2}/2$ . The results of figure 6 reveal that the vertical velocity of the bulge is smaller than that associated with free fall during most of its time evolution. Only during the last stages, when the volume accumulated inside the fluid bulge becomes larger than the volume of the liquid filaments connecting the bulge with the injector upstream and with the reservoir downstream, the free-fall law is approached. This fact suggests that the liquid bulge behaves like a freely falling liquid drop when its volume becomes large enough.

Figure 7. Bifurcation diagram for $\unicode[STIX]{x1D6E4}=16.67$ , $Bo=1$ and $L/R=23.26$ , where the square of the saturated amplitude of the radius oscillations, $A_{sat}^{2}$ (defined in the main text), is represented as a function of $We$ . Results from experiments (▵) and numerical integrations of (2.1)–(2.3) (○) are shown. The experimental critical Weber number lies in the range $5.62\times 10^{-3}<We_{c}^{exp}<5.83\times 10^{-3}$ , represented by the symbols ▸ and ◂. The numerical simulations provide the value $We_{c}^{num}=7.88\times 10^{-3}$ (●), in agreement with the result of the linear stability analysis, $We_{c}^{lsa}$ . The solid lines are linear fits close to the experimental and numerical neutral points, according to the supercritical Stuart–Landau model. The saturated oscillation frequencies, $\unicode[STIX]{x1D714}_{i}$ , are indicated close to several data points. The inset shows the structure of the numerical limit cycle in the $(A_{u},A_{r_{j}})$ plane associated with the grey-filled circle. The offset between the values of $We_{c}^{exp}$ and $We_{c}^{num}$ can also be seen in figures 2 and 9(d).

The bifurcation diagram represented in figure 7 shows the squared amplitude of the saturated radius oscillations, $A_{sat}^{2}$ , as a function of $We$ , obtained from the experiments (▵) and numerical simulations (○). Here the saturated amplitude is defined as $A_{sat}=\max (A_{r_{j}})-\min (A_{r_{j}})$ with $t>t_{sat}$ , $t_{sat}$ being a value of time large enough for the asymptotic limit-cycle behaviour illustrated by the inset of figure 7 to be reached. Moreover, the amplitude is defined by the same pointwise measure used in figure 3, namely $A_{r_{j}}=r_{j}(z^{\ast },t)-r_{j0}(z^{\ast })$ and, similarly, $A_{u}=u(z^{\ast },t)-u_{0}(z^{\ast })$ . The results shown in figure 7 prove that the axially confined jets under study behave as hydrodynamic oscillators governed by a supercritical Hopf bifurcation. Indeed, in the case of the numerical simulations, as the value of $We^{num}$ decreases and the critical value $We_{c}^{lsa}$ given by the linear analysis of figure 2 is crossed (●), a branch of finite-amplitude, orbitally stable periodic solutions emerges for $We^{num}<We_{c}^{lsa}$ (○). A similar scenario holds also for the experiments, but at smaller values of $We^{exp}$ , since $We_{c}^{exp}<We_{c}^{lsa}$ as deduced from figure 2. In addition, figure 7 shows that close to the critical point the amplitude grows as $A_{sat}^{2}\propto (We_{c}-We)$ , as expected from the Stuart–Landau equation (Landau Reference Landau1944; Stuart Reference Stuart1958). Note also that the frequency decreases as the value of $We$ decreases, whilst the growth rate increases leading to shorter saturation times.

Figure 8. (a) Photograph extracted from movie 3 of the supplementary material, illustrating the incipient breakup of a jet with $\unicode[STIX]{x1D708}=500~\text{cSt}$ , $R=1.75~\text{mm}$ , $L=16.6~\text{mm}$ , injected at a flow rate, $Q=3.3~\text{ml}~\text{min}^{-1}$ . These conditions are represented in the phase diagrams of figure 9(a,c) with the symbol ▴, which lies within the dripping region. Since $Q<Q_{b}$ , the oscillation amplitude increases and finally leads to dripping through the breakup of a thin thread that connects the meniscus attached to the injector with the liquid reservoir, as clearly observed in movie 3. (b) Instantaneous jet shape at incipient breakup obtained by numerically integrating equations (2.1)–(2.3) for the equivalent dimensionless configuration, namely $\unicode[STIX]{x1D6E4}=8.33$ , $Bo=1.38$ , $We=2.63\times 10^{-3}$ , $L/R=9.49$ . An animated numerical sequence under these conditions is displayed in movie 4 of the supplementary material.

Since, as revealed by figure 5, the saturated oscillation amplitude increases with $z$ , it can be anticipated that the breakup of the liquid thread will occur as the value of $L/R$ increases, leading to a jetting–dripping transition. Moreover, since the minimum jet radius is always located slightly upstream of the impact region, pinch-off should take place near the reservoir. Indeed, this fact is shown experimentally in the photograph of figure 8(a), extracted from movie 3 of the supplementary material, and numerically in figure 8(b), extracted from movie 4 of the supplementary material. Note that the one-dimensional model correctly captures the shape of the interface close to breakup except for the impact region downstream of the minimum radius, probably due to the crude model adopted herein of the outlet boundary condition.

Figure 9. Phase diagrams represented in the dimensional parameter plane $(Q,L)$ (a,b), and in the equivalent dimensionless parameter plane $(We,L/R)$ (c,d), showing the different regimes identified for an axially confined viscous liquid jet falling under gravity. The dark grey, light grey and unshaded regions correspond with the steady jetting, oscillatory jetting and dripping regimes, respectively. Results are shown for two different combinations of liquid and injector radius, namely (a,c) $\unicode[STIX]{x1D708}=500~\text{cSt}\;(\unicode[STIX]{x1D6E4}=8.33)$ and $R=1.75~\text{mm}$ ( $Bo=1.38$ ) and (b,d) $\unicode[STIX]{x1D708}=1000~\text{cSt}\;(\unicode[STIX]{x1D6E4}=16.67)$ , $R=1.5~\text{mm}\;(Bo=1)$ . The mean between the upper and lower limits of the experimental critical flow rate, $Q_{c}^{exp}$ ( $We_{c}^{exp}$ ), and breakup flow rate, $Q_{b}^{exp}$ ( $We_{b}^{exp}$ ), are plotted with the symbols ○ and ▫, respectively. The solid line is the critical flow rate for the onset of the breathing mode given by the linear stability analysis, $Q_{c}^{lsa}$ ( $We_{c}^{lsa}$ ), and the dashed line represents the function $Q_{b}^{num}$ ( $We_{b}^{num}$ ) obtained from the numerical integration of (2.1)–(2.3). The grey-filled circles correspond with direct transitions between steady jetting and dripping, while the symbols ▴ represents the conditions of figure 8. The different coiling states identified in the experiments are indicated.

The preceding discussion naturally leads to the following question: for a constant value of $We$ , how much can $L/R$ increase while avoiding the breakup of the liquid filament? Or, alternatively, how much can $We$ decrease for a constant value of $L/R$ to avoid pinch-off? To address these questions we will make use of the breakup flow rate $Q_{b}$ or, in dimensionless terms, the breakup Weber number, $We_{b}$ . The experimental and numerical phase diagrams represented in figure 9 show the long time regimes reached by the liquid jet in the $(Q,L)$ dimensional parameter plane in (a,b), and in the equivalent dimensionless parameter plane $(We,L/R)$ in (c,d). Two different combinations of liquid and injector radius are reported in figure 9, namely $\unicode[STIX]{x1D708}=500$ cSt ( $\unicode[STIX]{x1D6E4}=8.33$ ) and $R=1.75~\text{mm}\;(Bo=1.38)$ in (a,c), while $\unicode[STIX]{x1D708}=1000~\text{cSt}\;(\unicode[STIX]{x1D6E4}=16.67)$ and $R=1.5~\text{mm}\;(Bo=1)$ in (b,d). The dark and light grey regions correspond to the steady jetting and oscillatory jetting states according to the experiments, while the dripping regime falls within the unshaded regions. The solid curve is the critical flow rate $Q_{c}^{lsa}$ ( $We_{c}^{lsa}$ ) given by the linear stability analysis of § 3.1, and the dashed line is the breakup flow rate $Q_{b}^{num}$ ( $We_{b}^{num}$ ) according to the numerical simulations of (2.1)–(2.3), obtained by extrapolating the minimum saturated jet radius for several values of $We$ slightly larger than $We_{b}^{num}$ . The agreement between $We_{b}^{exp}$ and $We_{b}^{num}$ is fairly good for the case with $\unicode[STIX]{x1D6E4}=16.67\;(\unicode[STIX]{x1D708}=1000~\text{cSt})$ , but quite poor for the case with $\unicode[STIX]{x1D6E4}=8.33\;(\unicode[STIX]{x1D708}=500~\text{cSt})$ and values of $L/R\gtrsim 10$ . This discrepancy between theory and experiment for values of $L/R\gtrsim 10$ may be due to the limitations of the one-dimensional model that have already been pointed out in the discussion of figure 2. Indeed, the model does not account for the downstream viscous relaxation of the parabolic velocity profile at the injector outlet. More importantly, our crude modelling of the non-slender impact region where, as shown in figure 8, the breakup of the thread takes place, may also contribute to the difference between the values of $We_{b}^{exp}$ and $We_{b}^{num}$ observed in figure 9.

A salient feature of the phase maps is the existence of a critical length, $L^{\ast }(\unicode[STIX]{x1D708},R)$ for which $Q_{c}=Q_{b}=Q^{\ast }$ . For values of $L<L^{\ast }$ a direct transition between steady jetting and dripping takes place, without the existence of an intermediate oscillating jet state. The corresponding transition points are shown as grey-filled circles in figure 9, displaying a very good agreement with the value of $Q_{c}^{lsa}$ . It is important to emphasise that this phenomenon was observed both in the experiments and in the numerical simulations. Let us also point out that, in the cases where a direct transition from steady jetting to dripping takes place, there is no qualitative change in the linear stability mode, which is still the same breathing mode that is destabilised for $Q<Q_{c}$ . Moreover, the nature of the bifurcation that takes place in these cases is also the same Hopf bifurcation shown in figure 7, the only difference being that the amplitude of the oscillations grows with time without saturation, until breakup takes place, for very small values of $Q_{c}-Q$ . Ideally there should always exist a value of $Q_{c}-Q$ small enough for the oscillatory jetting regime to be achieved, since the amplitude of the oscillations increases as $(Q_{c}-Q)^{1/2}$ for sufficiently small values of $Q_{c}-Q$ . In practice, however, the value of $Q_{c}-Q_{b}$ is so small in these cases that it cannot be measured in our experiments or numerical simulations.

Figure 9 also reveals that the region of oscillatory jetting decreases as the value of $L/R$ increases, since the amplitude of the oscillations grows with $z$ , and thus $Q_{b}$ increases with $L$ , while the value of $Q_{c}$ reaches an asymptote as $L\rightarrow \infty$ (see also figure 2). Therefore, for values of $L$ larger than those considered herein, it is expected that the values of $Q_{c}$ and $Q_{b}$ will intersect, providing a bounded region of oscillatory jetting. Finally, figure 9 shows the important influence of $\unicode[STIX]{x1D708}$ , which increases the size of oscillatory jetting region due to its stabilising role. Similarly, the region of oscillatory jetting becomes smaller as $R$ increases due to the stabilising effect of axial stretching.