1 Introduction
The great research effort devoted in the past to understanding and controlling the dynamics of liquid jets is justified by their rich phenomenology and the large number of technological applications where they play a central role, such as fuel atomisation, chemical reactors, ink-jet and 3D printing, additive manufacturing, microfluidic platforms, drug encapsulation, mass spectrometry or cytometry, to name a few (see e.g. the reviews by Bogy Reference Bogy1979; Eggers Reference Eggers1997; Lin & Reitz Reference Lin and Reitz1998; Basaran Reference Basaran2002; Barrero & Loscertales Reference Barrero and Loscertales2007; Christopher & Anna Reference Christopher and Anna2007; Eggers & Villermaux Reference Eggers and Villermaux2008; Derby Reference Derby2010; Anna Reference Anna2016). These applications often require the generation of micrometre-sized jets and drops, in which case it proves convenient to downscale the disperse phase from a typically millimetre-sized injector to the desired micrometre scales to avoid clogging issues. This stretching effect can be achieved through different techniques like fibre spinning (Matovich & Pearson Reference Matovich and Pearson1969; Pearson & Matovich Reference Pearson and Matovich1969), electrospinning (Doshi & Reneker Reference Doshi and Reneker1995), flow focusing (Gañán-Calvo Reference Gañán-Calvo1998), viscous co-flows (Suryo & Basaran Reference Suryo and Basaran2006; Marín, Campo-Cortés & Gordillo Reference Marín, Campo-Cortés and Gordillo2009; Evangelio, Campo-Cortés & Gordillo Reference Evangelio, Campo-Cortés and Gordillo2016) and gravitational stretching (Sauter & Buggisch Reference Sauter and Buggisch2005; Rubio-Rubio, Sevilla & Gordillo Reference Rubio-Rubio, Sevilla and Gordillo2013). The latter configuration is particularly simple to implement, and the present work naturally extends the investigation of Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013) to assess the influence of a finite jet length and nonlinearity on the dynamics of the liquid thread.
In the absence of axial confinement, previous studies of the downwards injection of a constant flow rate of liquid into a passive gaseous atmosphere have demonstrated the existence of two different flow states, namely dripping and jetting, respectively prevailing for small and large enough values of the flow rate (Clanet & Lasheras Reference Clanet and Lasheras1999; Ambravaneswaran et al. Reference Ambravaneswaran, Subramani, Phillips and Basaran2004). The jetting regime is characterised by the formation of a liquid column which breaks up into drops at a certain distance from the injector due to the downstream growth of axisymmetric capillary instability waves (see e.g. Plateau Reference Plateau1873; Rayleigh Reference Rayleigh1878; Donnelly & Glaberson Reference Donnelly and Glaberson1966; Kalaaji et al. Reference Kalaaji, Lopez, Attane and Soucemarianadin2003; González & García Reference González and García2009). In contrast, the dripping regime features the emission of comparatively larger drops near the injector (Wilkes, Phillips & Basaran Reference Wilkes, Phillips and Basaran1999; Coullet, Mahadevan & Riera Reference Coullet, Mahadevan and Riera2005; Subramani et al. Reference Subramani, Yeoh, Suryo, Xu, Ambravaneswaran and Basaran2006). From the point of view of local stability theory, the jetting regime is a convectively unstable flow, in which the downstream advection of growing disturbances by the underlying base flow allows the formation of a long liquid column. Moreover, several works have successfully described the transition from jetting to dripping as a global instability that, in the case of quasi-parallel jets, is linked with the onset of local absolute instability (Leib & Goldstein Reference Leib and Goldstein1986a ,Reference Leib and Goldstein b ; Le Dizès Reference Le Dizès1997; Vihinen, Honohan & Lin Reference Vihinen, Honohan and Lin1997; O’Donnell, Chen & Lin Reference O’Donnell, Chen and Lin2001; Söderberg Reference Söderberg2003; Sevilla Reference Sevilla2011; Guerrero, González & García Reference Guerrero, González and García2016).
The concepts of convective and absolute instability rely on the assumption of quasi-parallel flow, and do not apply to cases where the wavelength of disturbances is of the order of the development length of the base flow, as happens in highly stretched jets like those studied in the present work. In these cases a global stability analysis must be performed, in which the spatial structure of the eigenfunctions is obtained as part of the solution (Theofilis Reference Theofilis2011). The global stability analysis is greatly simplified by the use of one-dimensional approximations to the full conservation equations, since the eigenvalue problem involves only the axial coordinate as an eigendirection (see e.g. Pearson & Matovich Reference Pearson and Matovich1969; Sauter & Buggisch Reference Sauter and Buggisch2005; Rubio-Rubio et al. Reference Rubio-Rubio, Sevilla and Gordillo2013; Gordillo, Sevilla & Campo-Cortés Reference Gordillo, Sevilla and Campo-Cortés2014). In particular, a global stability analysis of the leading-order one-dimensional model for viscous liquid columns (Eggers & Dupont Reference Eggers and Dupont1994; García & Castellanos Reference García and Castellanos1994) was first applied to gravitationally stretched viscous jets by Sauter & Buggisch (Reference Sauter and Buggisch2005), and refined later on by Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013). The latter works revealed that the axisymmetric self-excited oscillations observed in long viscous jets below a certain critical flow rate are explained by the destabilisation of a linear global mode.
Despite the usefulness of linearised theory to predict many features of liquid jets, there are relevant aspects of their dynamics that can only be described using a nonlinear approach, prominent examples being the pinch-off singularity (Eggers Reference Eggers1993) and the formation of satellite droplets (Yuen Reference Yuen1968; Nayfeh Reference Nayfeh1970; Rutland & Jameson Reference Rutland and Jameson1971; Ashgriz & Mashayek Reference Ashgriz and Mashayek1995). In the present work we introduce a new configuration where nonlinearity provides the selection mechanism between two different regimes after the jet becomes globally unstable: either a limit-cycle state without breakup described here for the first time (see movie 2 of the supplementary material available at https://doi.org/10.1017/jfm.2017.706), or a fully developed dripping state (see movie 3 of the supplementary material).
Liquid threads of finite length, either by their impingement onto a bath of the same liquid or by their impact on a solid plate, have also been considered in previous works. Thus, the shape and stability of quasi-static axisymmetric liquid bridges formed between a solid rod and an infinite bath were studied by Kovitz (Reference Kovitz1975), and more recently by Benilov & Oron (Reference Benilov and Oron2010) and Benilov & Cummins (Reference Benilov and Cummins2013). Christodoulides & Dias (Reference Christodoulides and Dias2010) studied the steady structure of planar vertical inviscid jets impinging onto a horizontal solid plate. However, most of the literature deals with the phenomenon of coiling, the fascinating buckling instability associated with the impact of a jet or film of viscous liquid on a solid substrate (see Ribe, Habibi & Bonn Reference Ribe, Habibi and Bonn2012, and references therein). Another beautiful and surprisingly rich phenomenon that has been studied is the deposition of viscous threads on moving solid substrates (Chiu-Webster & Lister Reference Chiu-Webster and Lister2006; Blount & Lister Reference Blount and Lister2011). These coiling states and deposition patterns break the axial symmetry, and their mathematical description is complicated by the fact that the thread centreline must be obtained as part of the solution (Entov & Yarin Reference Entov and Yarin1984). Although in the experiments reported herein we have observed coiling under certain conditions, our focus is on the global self-excited oscillations caused by the destabilisation of the axisymmetric breathing mode studied by Sauter & Buggisch (Reference Sauter and Buggisch2005) and Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013). It should be pointed out that the breathing mode and coiling coexist in a wide region of parameter space, as evidenced by movie 5 of the supplementary material. Nevertheless, in all the cases considered herein the coupling between both modes is weak enough for a purely axisymmetric model to provide a reasonably good leading-order description, not only of the neutral conditions for the onset of the breathing mode, but also of the long-time regime of the jet.
In the present work we report experiments performed to characterise the linear and nonlinear stability properties of jets of Newtonian liquid, injected at a constant flow rate through a circular tube, that impinge on the free surface of a reservoir of the same liquid placed at a controlled distance from the injector. The experiments are complemented with a global stability analysis and numerical simulations, both based on the leading-order one-dimensional mass and momentum equations retaining the full expression of the interfacial curvature (Eggers & Dupont Reference Eggers and Dupont1994; Rubio-Rubio et al. Reference Rubio-Rubio, Sevilla and Gordillo2013). We also study the selection of nonlinear regimes under globally unstable conditions. The experimental set-up and the mathematical model are described in § 2, and the results are presented in § 3, finally leading to the conclusions drawn in § 4.
2 Experimental set-up and mathematical model
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.
(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.
-
(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}$
:
-
(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.
-
(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),
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}r_{j}^{2}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}(ur_{j}^{2})}{\unicode[STIX]{x2202}z}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}t}+u\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}=1-\frac{\unicode[STIX]{x2202}{\mathcal{C}}}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x1D6E4}}{r_{j}^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left(r_{j}^{2}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{C}}=r_{j}^{-1}\left[1+\left(\frac{\unicode[STIX]{x2202}r_{j}}{\unicode[STIX]{x2202}z}\right)^{2}\right]^{-1/2}-\frac{\unicode[STIX]{x2202}^{2}r_{j}}{\unicode[STIX]{x2202}z^{2}}\left[1+\left(\frac{\unicode[STIX]{x2202}r_{j}}{\unicode[STIX]{x2202}z}\right)^{2}\right]^{-3/2}, & \displaystyle\end{eqnarray}$$
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 Experimental and theoretical results
In this section we present the results of the experiments, as well as the linear stability analysis and the numerical simulations, based on the one-dimensional model (2.1)–(2.3).
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,
$$\begin{eqnarray}\displaystyle & \displaystyle r_{j}(z,t)=r_{j0}(z)+\unicode[STIX]{x1D716}r_{j1}(z)\text{e}^{\unicode[STIX]{x1D714}t}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u(z,t)=u_{0}(z)+\unicode[STIX]{x1D716}u_{1}(z)\text{e}^{\unicode[STIX]{x1D714}t}, & \displaystyle\end{eqnarray}$$
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).
(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.
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$
.
(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.
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.
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.
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.
(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.
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.
4 Conclusions
The linear and nonlinear dynamics of jets of viscous liquid falling under gravity and confined in the axial direction has been studied by means of experiments and a simple one-dimensional model (Eggers & Dupont Reference Eggers and Dupont1994; García & Castellanos Reference García and Castellanos1994) that is extremely efficient from a computational point of view.
The global linear stability analysis reported by Sauter & Buggisch (Reference Sauter and Buggisch2005) and Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013) has been extended in the present work to contemplate the effect of the axial confinement length,
$L$
. We have shown both experimentally and theoretically that the critical flow rate
$Q_{c}$
for the onset of the axisymmetric breathing mode decreases as
$L$
decreases, thus stabilising the liquid thread. For small enough values of
$L$
, a limit is reached whereby
$We_{c}\rightarrow 0$
and the marginal stability point is given by hydrostatics, in agreement with previous works on static liquid bridges (Kovitz Reference Kovitz1975; Benilov & Oron Reference Benilov and Oron2010; Benilov & Cummins Reference Benilov and Cummins2013).
The nonlinear dynamics has been characterised both experimentally and with numerical integrations of the one-dimensional model. Our results have revealed the existence of a new regime featuring the appearance limit-cycle oscillations of the thread without breakup. The latter oscillatory jetting regime has been thoroughly studied for several combinations of the governing parameters, finally leading to the phase maps presented in figure 9, that summarises the main conclusions of the present work. In particular, the central role of
$L$
is clearly unveiled. In contrast with the linear theory, which only distinguishes between stable and unstable configurations, three distinct long-time nonlinear regimes can be observed in figure 9: steady jetting, oscillatory jetting and dripping. Although different coiling regimes have also been identified in the experiments, namely a standard and an oscillatory coiling, a detailed study of these coiling regimes is out of the scope of the present work.
The experiments and numerical simulations have revealed that there is a critical value of the confinement length,
$L$
, below which the oscillatory jetting is no longer observed. This is due to the fact that, as
$L$
decreases,
$Q_{c}$
decreases with a larger slope than the critical flow rate associated with the breakup of the jet,
$Q_{b}$
. Hence, the size of the oscillatory jetting region decreases as the jet becomes shorter, due to the stabilising effect of confinement. Consequently, there is a value of
$L=L^{\ast }$
where
$Q_{c}=Q_{b}=Q^{\ast }$
, and for values of
$L<L^{\ast }$
the transition occurs directly between steady jetting and dripping as
$Q$
decreases. The experiments and numerical simulations are in very good agreement in their prediction of the value of
$L^{\ast }$
and the corresponding
$Q^{\ast }$
.
Future work should contemplate wider ranges of the control parameters, allowing to obtain experimental and numerical phase maps for values of
$\unicode[STIX]{x1D6E4}$
and
$Bo$
different from those reported herein. Additionally, in view of the new coiling regimes found in our experiments, it could be worthwhile to describe them in detail. Another important extension is the characterisation of the physical effects present in typical technological applications, such as the use of liquids with complex rheology or the presence of surfactants. Their inclusion in the analysis is crucial for applications such as 3D printing or additive manufacturing, in which our findings may be of relevance.
Acknowledgements
The authors thank the Spanish MINECO, Subdirección General de Gestión de Ayudas a la Investigación, for its support through projects DPI2014-59292-C3-1-P, DPI2014-59292-C3-3-P and DPI2015-71901-REDT. This research project has been partly financed through European funds.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2017.706.
Appendix A. Linear stability analysis
The detailed mathematical formulation of the linear stability problem is provided in the present appendix.
A.1 Base flow
The base flow [
$r_{j0}(z)$
,
$u_{0}(z)$
] satisfies the steady versions of the system of (2.1)–(2.3). In particular, the continuity equation (2.1) allows us to substitute the steady state velocity,
$u_{0}(z)$
, by the base flow jet radius,
$r_{j0}(z)$
,
$$\begin{eqnarray}u_{0}=\frac{q}{r_{j0}^{2}},\end{eqnarray}$$
where
$q=We^{1/2}Bo^{3/4}$
, is the dimensionless liquid flow rate. In addition, the steady version of the momentum equation (2.2) reduces to
$$\begin{eqnarray}u_{0}{u_{0}}^{\prime }=1-{{\mathcal{C}}_{0}}^{\prime }+\unicode[STIX]{x1D6E4}{r_{j0}}^{-2}({r_{j0}}^{2}{u_{0}}^{\prime })^{\prime },\end{eqnarray}$$
where primes indicate derivatives with respect to
$z$
and
${\mathcal{C}}_{0}$
is the interfacial curvature associated with the steady jet shape,
$r_{j0}(z)$
, which satisfies equation (2.3). Substituting (A 1) into (A 2), the function
$r_{j0}(z)$
satisfies
$$\begin{eqnarray}-r_{j0}^{2}{{\mathcal{C}}_{0}}^{\prime }+2q\unicode[STIX]{x1D6E4}[r_{j0}^{-2}(r_{j0}^{\prime })^{2}-r_{j0}^{-1}r_{j0}^{\prime \prime }]+2q^{2}r_{j0}^{-3}r_{j0}^{\prime }+r_{j0}^{2}=0,\end{eqnarray}$$
where
$$\begin{eqnarray}-r_{j0}^{2}{{\mathcal{C}}_{0}}^{\prime }=\frac{r_{j0}^{\prime }}{[1+(r_{j0}^{\prime })^{2}]^{1/2}}+\frac{r_{j0}r_{j0}^{\prime }r_{j0}^{\prime \prime }+r_{j0}^{2}r_{j0}^{\prime \prime \prime }}{[1+(r_{j0}^{\prime })^{2}]^{3/2}}-\frac{3r_{j0}^{2}r_{j0}^{\prime }(r_{j0}^{\prime \prime })^{2}}{[1+(r_{j0}^{\prime })^{2}]^{5/2}},\end{eqnarray}$$
to be solved with the boundary conditions
$$\begin{eqnarray}\displaystyle & \displaystyle z=0:r_{j0}=Bo^{1/2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle z=0:u_{0}=We^{1/2}Bo^{-1/4}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle z=\frac{L}{l_{\unicode[STIX]{x1D70E}}}:u_{0}=u_{out}, & \displaystyle\end{eqnarray}$$
as discussed in § 2.3. Note that (A 3)–(A 4) are the same as those solved in Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013) (hereinafter R13), the only difference being the outlet boundary condition (A 7), which in the present work is imposed as a Dirichlet condition, while a free boundary condition was considered in R13. The numerical method used herein to solve equations (A 3)–(A 4) is also the same as that used by R13, and is explained in § B.1.
A.2 Linear stability problem
Substituting (3.1)–(3.2) into (2.1)–(2.3), the
$O(\unicode[STIX]{x1D716})$
eigenvalue problem can be written in the following compact form,
$$\begin{eqnarray}\left[\begin{array}{@{}cc@{}}{\mathcal{M}}_{r_{j}}^{c} & {\mathcal{M}}_{u}^{c}\\ {\mathcal{M}}_{r_{j}}^{m} & {\mathcal{M}}_{u}^{m}\end{array}\right]\left[\begin{array}{@{}c@{}}r_{j1}\\ u_{1}\end{array}\right]=\unicode[STIX]{x1D714}\left[\begin{array}{@{}c@{}}r_{j1}\\ u_{1}\end{array}\right],\end{eqnarray}$$
with
${\mathcal{M}}_{i}^{j}$
denoting the following differential operators,
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{M}}_{r_{j}}^{c}=-\frac{q}{r_{j0}^{2}}D+\frac{qr_{j0}^{\prime }}{r_{j0}^{3}}I, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{M}}_{u}^{c}=-\frac{r_{j0}}{2}D-r_{j0}^{\prime }I, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{M}}_{r_{j}}^{m}=-4\unicode[STIX]{x1D6E4}q\left(\frac{r_{j0}^{\prime }}{r_{j0}^{4}}D-\frac{(r_{j0}^{\prime })^{2}}{r_{j0}^{5}}I\right)+\mathop{\sum }_{k=1}^{4}{\mathcal{S}}^{2k-1}{\mathcal{H}}_{k}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{M}}_{u}^{m}=\unicode[STIX]{x1D6E4}\left(D^{2}+\frac{2r_{j0}^{\prime }}{r_{j0}}D\right)-\frac{q}{r_{j0}^{2}}D+\frac{2qr_{j0}^{\prime }}{r_{j0}^{3}}I. & \displaystyle\end{eqnarray}$$
In (A 9)–(A 12),
$I$
is the identity,
$D^{n}\equiv \text{d}^{n}/\text{d}z^{n}$
,
${\mathcal{S}}(z)=[1+(r_{j0}^{\prime })^{2}]^{-1/2}$
and
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{H}}_{1}=\frac{1}{r_{j0}^{2}}D-\frac{2r_{j0}^{\prime }}{r_{j0}^{3}}I, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{H}}_{2}=D^{3}+\frac{r_{j0}^{\prime }}{r_{j0}}D^{2}-\left[\frac{(r_{j0}^{\prime })^{2}}{r_{j0}^{2}}-\frac{r_{j0}^{\prime \prime }}{r_{j0}}\right]D-\frac{r_{j0}^{\prime }r_{j0}^{\prime \prime }}{r_{j0}^{2}}I, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{H}}_{3}=-6r_{j0}^{\prime }r_{j0}^{\prime \prime }D^{2}-3\left[\frac{(r_{j0}^{\prime })^{2}r_{j0}^{\prime \prime }}{r_{j0}}+(r_{j0}^{\prime \prime })^{2}+r_{j0}^{\prime }r_{j0}^{\prime \prime \prime }\right]D, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{H}}_{4}=15(r_{j0}^{\prime \prime })^{2}(r_{j0}^{\prime })^{2}D. & \displaystyle\end{eqnarray}$$
Note that (A 8) is a linear eigenvalue problem, complemented with the boundary conditions
$$\begin{eqnarray}\displaystyle & \displaystyle z=0:u_{1}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle z=0:r_{j1}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle z=\frac{L}{l_{\unicode[STIX]{x1D70E}}}:u_{1}=0. & \displaystyle\end{eqnarray}$$
Equations (A 17) and (A 18) represent the pinned contact line and constant flow rate conditions, respectively. However, equation (A 19) just represents in a crude way the fact that the disturbed jet velocity is small in the impact region onto the bath at rest. If a small enough value of
$u_{out}$
is assumed for the base flow, the results of the stability analysis do not vary significantly, thereby justifying the use of (A 19). Moreover, for values of the jet length
$L/R\gtrsim 20$
the results are independent of the value of
$u_{out}$
, and in fact we have checked that both the base flow and its linear stability are virtually the same whether imposing a Neumann condition at
$z=L/l_{\unicode[STIX]{x1D70E}}$
, or even not imposing any condition at all, in agreement with the results of R13.
Appendix B. Numerical methods
In the present appendix we describe the numerical methods used for the computation of the steady jet and its linear stability, as well as the integration of (2.1)–(2.3).
B.1 Base flow and linear stability analysis
To solve both (A 3) for
$r_{j0}$
, and the system (A 8) for the eigenvalues
$\unicode[STIX]{x1D714}$
and the corresponding eigenfunctions, we discretised the differential operators using a Chebyshev collocation method (Canuto et al.
Reference Canuto, Hussaini, Quarteroni and Zang2006). To that end, the physical domain,
$0\leqslant z\leqslant L/l_{\unicode[STIX]{x1D70E}}$
is mapped into the interval
$-1\leqslant y\leqslant 1$
by means of the transformation
$$\begin{eqnarray}z=\frac{bL/l_{\unicode[STIX]{x1D70E}}(1+y)}{2b+L/l_{\unicode[STIX]{x1D70E}}(1-y)},\end{eqnarray}$$
where
$b$
is a parameter that controls the clustering of nodes at
$z=0$
and
$z=L/l_{\unicode[STIX]{x1D70E}}$
. Derivatives with respect to
$z$
are calculated using the standard Chebyshev differentiation matrices and the chain rule. The nonlinear differential equation (A 3) with boundary conditions (A 5)–(A 7) is solved first using an iterative Newton–Raphson method. Once
$r_{j0}$
is known at the
$N$
Chebyshev collocation points, the discretised version of (A 8), which results in a linear algebraic eigenvalue problem to determine the
$2N$
eigenvalues
$\unicode[STIX]{x1D714}^{k}$
,
$k=1\cdots 2N$
, and their corresponding eigenfunctions
$(r_{j1}^{k},u_{1}^{k})$
at the
$N$
collocation points, is solved using standard Matlab routines. Notice that the eigenfunctions must satisfy the boundary conditions (A 17)–(A 19). Although all the results reported in the present paper were computed with values of
$N$
and
$b$
within the ranges
$30\leqslant N\leqslant 200$
and
$5\leqslant b\leqslant 80$
, we have carefully checked that the leading eigenvalues and eigenfunctions are insensitive to the values of these parameters.
B.2 Direct numerical simulation
The system of (2.1)–(2.3), supplemented with the boundary conditions discussed above were solved with a very simple and efficient method of lines in which the spatial derivatives were approximated using the Chebyshev collocation method described in § B.1. The resulting system of
$2N$
nonlinear ordinary differential equations were integrated in time with an initial condition taken as the base flow
$[r_{j0}(z),u_{0}(z)]$
, slightly perturbed by a Gaussian function of very small amplitude. The ode23t routine from the Matlab ODE suite was chosen, since Dirichlet or Neumann boundary conditions can be easily imposed through a mass matrix, and also allows us to implement the Jacobian matrix of the nonlinear system to improve the speed and accuracy of the computations. The solver ode23t uses the Bogacki–Shampine algorithm, which is a method based on two single-step formulas, of second and third order respectively, and computes three stages for each time step. Briefly stated, to calculate the temporal evolution of the
$r_{j}$
and
$u$
, the semi-discretised equations were written in the matrix form
$$\begin{eqnarray}\frac{\text{d}\boldsymbol{y}}{\text{d}t}=\unicode[STIX]{x1D63D}(\boldsymbol{y},t),\end{eqnarray}$$
where the column vector
$\boldsymbol{y}=[\boldsymbol{r}_{j}(\boldsymbol{z},t),\boldsymbol{u}(\boldsymbol{z},t)]^{\text{T}}$
contains the values of the jet radius and axial velocity computed at the
$N$
Chebyshev collocation points,
$\boldsymbol{z}$
, and the matrix
$\unicode[STIX]{x1D63D}$
is
The
$N\times N$
submatrices
$\unicode[STIX]{x1D63D}_{i}^{j}$
appearing in (B 3) are
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63D}_{r_{j}}^{c}=-\text{diag}(\boldsymbol{u})\boldsymbol{\cdot }\unicode[STIX]{x1D63F}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63D}_{u}^{c}=-\frac{1}{2}\,\text{diag}(\boldsymbol{r}_{j})\boldsymbol{\cdot }\unicode[STIX]{x1D63F}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63D}_{r_{j}}^{m}=\mathop{\sum }_{n=1}^{3}\unicode[STIX]{x1D63E}_{n}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}^{n}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63D}_{u}^{m}=-\text{diag}(\boldsymbol{u})\boldsymbol{\cdot }\unicode[STIX]{x1D63F}+\unicode[STIX]{x1D6E4}[\unicode[STIX]{x1D63F}^{2}+2\,\text{diag}(\boldsymbol{r}_{j})^{-1}\boldsymbol{\cdot }\text{diag}(\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\boldsymbol{r}_{j})\boldsymbol{\cdot }\unicode[STIX]{x1D63F}], & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D63F}^{n}$
is the
$N\times N~n$
th order Chebyshev differentiation matrix,
$\text{diag}(\cdot )$
maps an
$N$
-tuple to the corresponding diagonal matrix, and the curvature gradient diagonal matrices,
$\unicode[STIX]{x1D63E}_{n}$
, are
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63E}_{1}=\text{diag}(\boldsymbol{r}_{j})^{-2}\boldsymbol{\cdot }[\unicode[STIX]{x1D644}+\text{diag}(\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\boldsymbol{r}_{j})^{2}]^{-1/2}-3\,\text{diag}(\unicode[STIX]{x1D63F}^{2}\boldsymbol{\cdot }\boldsymbol{r}_{j})^{2}\boldsymbol{\cdot }[\unicode[STIX]{x1D644}+\text{diag}(\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\boldsymbol{r}_{j})^{2}]^{-5/2},\quad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63E}_{2}=\text{diag}(\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\boldsymbol{r}_{j})\boldsymbol{\cdot }\text{diag}(\boldsymbol{r}_{j})^{-1}\boldsymbol{\cdot }[\unicode[STIX]{x1D644}+\text{diag}(\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\boldsymbol{r}_{j})^{2}]^{-3/2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63E}_{3}=[\unicode[STIX]{x1D644}+\text{diag}(\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\boldsymbol{r}_{j})^{2}]^{-3/2}, & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D644}$
is the
$N\times N$
identity matrix. Note that, for clarity, in (B 4)–(B 10) the functional dependence of
$\boldsymbol{r}_{j}$
and
$\boldsymbol{u}$
on
$(\boldsymbol{z},t)$
has been omitted. The boundary conditions (A 5)–(A 7) described in § 2.3 can be readily implemented by imposing
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{r}_{j}}{\text{d}t}(\boldsymbol{z}=0,t)=\frac{\text{d}\boldsymbol{u}}{\text{d}t}(\boldsymbol{z}=0,t)=\frac{\text{d}\boldsymbol{u}}{\text{d}t}(\boldsymbol{z}=L/l_{\unicode[STIX]{x1D70E}},t)=0. & \displaystyle\end{eqnarray}$$
