Frequency selection in a gravitationally stretched capillary jet in the jetting regime

TL;DR

This study investigates how external periodic forcing affects a gravitationally stretched capillary jet, revealing that the optimal forcing frequency depends on amplitude and can be accurately predicted by resolvent analysis.

Contribution

It introduces a combined use of nonlinear simulations and resolvent analysis to predict optimal forcing frequencies in a stretching capillary jet, improving understanding of jet breakup dynamics.

Findings

Optimal forcing frequency varies with forcing amplitude.

Resolvent analysis accurately predicts jet breakup length.

WKBJ formalism fails to match nonlinear simulation results.

Abstract

A capillary jet falling under the effect of gravity continuously stretches while thinning downstream. We report here the effect of external periodic forcing on such a spatially varying jet in the jetting regime. Surprisingly, the optimal forcing frequency producing the most unstable jet is found to be highly dependent on the forcing amplitude. Taking benefit of the one-dimensional Eggers & Dupont (J. Fluid Mech., vol. 262, 1994, 205-221) equations, we investigate the case through nonlinear simulations and linear stability analysis. In the local framework the WKBJ formalism, established for weakly non-parallel flows, fails to capture the nonlinear simulation results quantitatively. However in the global framework, the resolvent analysis supplemented by a simple approximation of the required response norm inducing breakup, is shown to correctly predict the optimal forcing frequency at a given forcing amplitude and the resulting jet breakup length. The results of the resolvent analysis are found to be in good agreement with those of the nonlinear simulations.

Equations

$$\frac{\partial h}{\partial t}$$

$$\frac{\partial u}{\partial t}$$

$$\phantom{-}p=\left(\frac{1}{h\Big{[}1+\Big{(}\frac{\partial h}{\partial z}\Big{)}^{2}\Big{]}^{\frac{1}{2}}}-\frac{\frac{\partial^{2}h}{\partial z^{2}}}{\Big{[}1+\Big{(}\frac{\partial h}{\partial z}\Big{)}^{2}\Big{]}^{\frac{3}{2}}}\right).$$

$$h(0,t)$$

$$u(0,t)$$

$$h_b^2{u}_b=Q=\sqrt{W{e}_{in}}.$$

$$2Q^2{h}_b^{\prime }-h_b^5{C}^{\prime }+6QO{h}_{in}(h_b{h_b^{\prime }}^2-h_b^2{h}_b^{\prime \prime })+h_b^{5}B{o}_{in}=0,$$

$$-h_b^2{C}^{\prime }=\frac{h_b^{\prime }}{[1+(h_b^{\prime }{)}^2{]}^{1/2}}+\frac{h_b{h}_b^{\prime }{h}_b^{\prime \prime }+h_b^2{h}_b^{\prime \prime \prime }}{[1+(h_b^{\prime }{)}^2{]}^{3/2}}-\frac{3h_b^2{h}_b^{\prime }(h_b^{\prime \prime }{)}^2}{[1+(h_b^{\prime }{)}^2{]}^{5/2}}.$$

$$\frac{\partial a}{\partial t}$$

$$\frac{\partial u}{\partial t}$$

$$p$$

$$u_f(0,t)=\mathrm{I}\mathrm{R}(ϵe^{\text{i}\omega t}),$$

$${\left(\frac{dA}{dz}\right)}_i=\beta {\left(\frac{dA}{dz}\right)}_{i,b}+(1-\beta ){\left(\frac{dA}{dz}\right)}_{i,f},$$

$$\beta =\frac{\tanh (\alpha u_i)+1}{2}.$$

$$h(z,t)$$

$$u(z,t)$$

$$a(z,t)=a_b+ϵ\hat{a}{e}^{\mathrm{i}(kz-\omega t)},$$

$$\omega =\sqrt{We}k-\frac{3\mathrm{i}Oh{k}^2}{2}\pm \mathrm{i}\sqrt{\frac{(k^2-k^4)}{2}+\frac{9{Oh}^2{k}^4}{4}},$$

$$\frac{d\Delta }{dk}({\omega }_0,k_0)=0,\Delta ({\omega }_0,k_0)=0,$$

$$\frac{\partial \Delta }{\partial k}=\sqrt{We}-3\mathrm{i}Ohk\pm \mathrm{i}\frac{1-2k^2+(3Ohk{)}^2}{\sqrt{2(1-k^2)+(3Ohk{)}^2}}.$$

$$O{h}_z$$

$$W{e}_z$$

$$\omega^{2}-2u_{b}(z)\omega k+\Bigg{(}\frac{1}{2\sqrt{a_{b}(z)}}+u_{b}(z)^{2}+3\text{i}\mathit{Oh_{z}}\omega\Bigg{)}k^{2}-3\text{i}\mathit{Oh_{z}}u_{b}(z)k^{3}-\frac{\sqrt{a_{b}(z)}}{2}k^{4}=0.$$

$$\mathbf{q^{\prime}}(z,t)=\mathbf{\hat{q}}(\omega,z)\exp\Bigg{[}\text{i}\Bigg{(}\int_{0}^{z}k(\omega,z^{\prime})dz^{\prime}-\omega t\Bigg{)}\Bigg{]}.$$

$$\mathit{G}_{s}(\mathit{\omega,L})=\left\lVert\displaystyle\int_{0}^{L}\mathbf{\hat{q}}(\omega,z)\exp\Bigg{[}\text{i}\Bigg{(}\int_{0}^{z}k(\omega,z^{\prime})dz^{\prime}-\omega t\Bigg{)}\Bigg{]}\right\rVert.$$

$$G_{s,max}(L)=\underset{\omega }{\max }[G_s(\omega ,L)],$$

$$\mathbf{q^{\prime}}(Z,t)\sim A(Z)\mathbf{\hat{q}}(Z,0)\exp\Bigg{[}\text{i}\Bigg{(}\frac{1}{\eta}\int_{0}^{Z}k(\omega,Z^{\prime})dZ^{\prime}-\omega t\Bigg{)}\Bigg{]},$$

$$M(Z)\frac{\mathrm{d}A(Z)}{\mathrm{d}Z}+N(Z)A(Z)=0,$$

$$A(Z)=A_{0}\exp\Bigg{(}-\int_{0}^{Z}\frac{N(Z^{\prime})}{M(Z^{\prime})}\mathrm{d}Z^{\prime}\Bigg{)}.$$

$$\mathit{G}_{Amp}^{2}(\mathit{\omega,L})=\frac{\displaystyle\int_{0}^{z}A^{H}(z^{\prime})A(z^{\prime})\Big{(}\mathbf{\hat{q}}^{H}(z^{\prime})\cdot\mathbf{\hat{q}}(z^{\prime})\Big{)}\Big{(}e^{\int_{0}^{z^{\prime}}-2k_{i}(z^{\prime\prime})\mathrm{d}z^{\prime\prime}}\Big{)}\mathrm{d}z^{\prime}}{\mathbf{\hat{q}}^{H}(0)\cdot\mathbf{\hat{q}}(0)}.$$

$$h(z,t)$$

$$u(z,t)$$

$$\lambda \mathrm{I}\mathrm{I}\left[\begin{array}{c}\tilde{h}\\ \tilde{u}\end{array}\right]=M\left[\begin{array}{c}\tilde{h}\\ \tilde{u}\end{array}\right],$$

$$\mathrm{I}\mathrm{I}{\partial }_{\mathrm{t}}[\mathbf{s}]=M[\mathbf{s}]+B_{f}f,$$

Full text

Physically, the breakup of the jet into drops begins with the excitation of a temporally or spatially amplifying suitable mode due to weak external disturbances. In practice, this weak agitation is usually imposed by controlled harmonic perturbation, either from within or at the outlet of the nozzle, to generate spatially amplifying waves leading to jet breakup. In this direction, the primary objective of this paper is to evaluate the response of an incompressible jet falling in presence of gravity, to externally imposed harmonic perturbations characterised by a fixed frequency and amplitude, and to find the optimal forcing which generates the most amplified response.

The external forcing is vital in the production of controlled micron-sized droplets, a feature essential to several application as in inkjet printers (Basaran 2002; Wijshoff 2010; Basaran et al. 2013), pharmaceuticals (Bennett et al., 2002) and powder technology (van Deventer et al., 2013), to name a few. In view of the limitations linked to the fabrication of such small droplets, most of the devices used for the drops production depend on the generation of highly thin liquid threads whose diameters are several orders smaller than the nozzle diameter. Some common methods for producing such threads use tangential electrical stresses (in electrospinning devices, Doshi & Reneker 1995; Loscertales et al. 2002), using outer co-flows (Marín et al., 2009) or a rotating spinerett (in fibre spinning applications (Pearson & Matovich, 1969)). Rubio-Rubio et al. (2013) showed an alternative method for producing highly elongated jets through the use of gravity, in which the mass conservation of the liquid jet forces its thinning as the liquid accelerates downstream.

The breakup of a liquid thread into drops, governed by the relative strength of the surface tension effect over the viscous and the inertial effects, was first explained by Plateau (1873) and Rayleigh (1879) for a uniform column of fluid. What adds complexity to the well understood viscous jet breakup mechanism is the presence of gravity which significantly stretches the base flow shape. The stability of the such spatially varying gravity jets should ideally be examined using the global stability analysis and by including the non-parallel effects of the base flow. A similar difficulty linked to the non-parallel nature of the flow results from the adaptation of the flow from a wall bounded flow within the nozzle to a free jet (Sevilla, 2011). Turning back to falling jets stretched by gravity, Sauter & Buggisch (2005) were the first to approach the problem theoretically by defining a linear global mode that correlated with the self sustained oscillations of the falling jet, observed during the jetting (globally stable) to dripping (globally unstable) transition. The work of Sauter & Buggisch (2005) was extended by Rubio-Rubio et al. (2013) experimentally and theoretically by increasing the range of liquid viscosities and nozzle diameters. Additionally they retained the entire expression of the curvature term for the formulation of their stability analysis, a feature that helped them to accurately predict the critical flow rate for the stability transition and the oscillating mode compared to the previous author. However, none of these studies predicted the jet stable length as a function of the flow rate and fluid properties, a question which was pursued by Javadi et al. (2013) experimentally and theoretically.

More recently, Le Dizès & Villermaux (2017) determined theoretically the stable jet length, wavelength at breakup and resulting drop size due to the most dangerous perturbation applied either at nozzle exit or affecting the jet all along its length for different jet viscosities. Their analysis accounted for both the base state deformation and modification of local instability dispersion relation as the jet thins in the direction of gravity. Notably, extending the work of previous authors (Tomotika, 1936; Frankel & Weihs, 1985; Leib & Goldstein, 1986; Frankel & Weihs, 1987; Senchenko & Bohr, 2005; Sauter & Buggisch, 2005; Javadi et al., 2013) they used the local plane wave decomposition (WKBJ approximation) for their analysis. However, the gain resulting from a perturbation was computed by considering only the exponential (e) terms of the WKBJ approximation. Additionally, an ad hoc spatial gain of $\text{e}^{7}$, of the linear perturbations was assumed to be sufficient for breakup. Thus the level of noise was considered fixed for all the theoretical analysis.

In this paper, we go beyond the global stability analysis of the gravity jets, and always operate in the stable regime where the jet behaves inherently as an amplifier. Precisely, we look at the receptivity of the jet in this regime to external perturbations, through nonlinear simulations and resolvent analysis with the aim of finding the optimal forcing which results in the most amplified disturbance. Unlike Le Dizès & Villermaux (2017), we consider an external forcing characterised by different amplitudes. The effect of forcing amplitude on the breakup length of very high speed jets has been numerically analysed by Hilbing & Heister (1996). However, a clear understanding of its effect in the case of spatially varying jets is still missing. Our analysis exemplifies the effect of forcing amplitude on the breakup length and the optimal forcing frequency. We also investigate the jet response using the WKBJ approximation and assess its validity for the spatially varying gravity jet. Our entire study is based on the slender-jet approximation (Eggers & Dupont, 1994) of the Navier-Stokes equation for an axisymmetric jet. The reduced one dimensional (1D) model has turned out to be extremely valuable for realistic representation of jets (Ambravaneswaran et al. 2002; van Hoeve et al. 2010) by accurately capturing the jet interface close to the breakup as well as the formation of ‘satellite’ drops. During the final stage of this work, we became aware of the work of (Lizzi, 2016) who has compared a resolvent analysis to experimental results and also revealed the dependance of the optimal frequency and breakup length as a function of time-harmonic forcing or noise amplitude.

The paper is structured as follows. Section 1 describes the governing equations. §2 discusses the nonlinear simulations where the results are detailed in §2.3. The local stability analysis of the gravity jet is performed in §3 where we compare the jet response using stability analysis in §3.3 and the WKBJ formulation in §3.4. We then operate in the global framework in §4 where the significance of the resolvent analysis is elucidated in §4.2. We show that the resolvent analysis is self-sufficient in predicting the optimal forcing frequency and the breakup length as obtained through the nonlinear simulations. Finally, we apply a white noise disturbance on the jet inlet to explore its behaviour in comparison to the expected response to the optimal forcing in §5. The conclusion and some perspectives related to the present work are summarised in §6.

1 Mathematical formulation

We consider an axisymmetric viscous jet falling vertically from a nozzle under the effect of gravity $g$. At the nozzle outlet, the jet has a fixed radius $\bar{h}_{0}$ and velocity $\bar{u}_{0}$. The surrounding medium is considered evanescent and is neglected. The density, dynamic viscosity and surface tension of the jet are denoted by $\rho$, $\mu$ and $\gamma$, respectively.

The behaviour of the jet is analysed using the leading-order one-dimensional mass and momentum equations, derived by Eggers & Dupont (1994). The dimensionless form of the equations, obtained by chosing $\bar{h}_{0}$ as the characteristic length scale, the inertial time $\tau_{i}=\sqrt{\rho\bar{h}_{0}^{3}/\gamma}$ as the characteristic time scale and $\gamma/\bar{h}^{2}_{0}$ as the pressure scale, are written as,

[TABLE]

where, the dimensionless pressure $p(z,t)$ is expressed as,

[TABLE]

In equation (1), $h(z,t)$ and $u(z,t)$ represent the height of the jet interface and the velocity at the axial distance $z$. The system of equations (1) are governed by the dimensionless numbers Ohnesorge ($\mathit{Oh_{in}}$) and Bond ($\mathit{Bo_{in}}$) defined at the inlet. Ohnesorge, expressed as \mathit{Oh_{in}}={\mu}\big{/}{\sqrt{\rho\gamma\bar{h}_{0}}}, relates the viscous forces to inertial and surface tension forces. The Bond (Eötvös) number denoted by \mathit{Bo_{in}}=\rho g{\bar{h}_{0}}^{2}\big{/}\gamma, measures the strength of the surface tension forces to body forces. A high $\mathit{Oh_{in}}$ or $\mathit{Bo_{in}}$ leads to a stabilised jet interface thus increasing the stability of the base flow.

Using the associated characteristic velocity $\bar{h}_{0}/\tau_{i}$, the non-dimensional boundary conditions for the jet at nozzle inlet are reduced to:

[TABLE]

Here, $\mathit{We_{in}}$ represents the Weber number defined at the nozzle inlet, \mathit{We_{in}}=\rho\bar{h}_{0}\bar{u}_{0}^{2}\big{/}\gamma, and measures the ratio between the kinetic energy and the surface energy.

The steady state form of the continuity equation (1a) gives the relation between the steady state shape $h_{b}$ and velocity $u_{b}$ as,

[TABLE]

Where $Q$ is the dimensionless flow rate, obtained from the nozzle conditions. This gives $u_{b}=\sqrt{\mathit{We_{in}}}/h_{b}^{2}$. Using the relation (4), the steady state momentum equation (1b) reduces to,

[TABLE]

where derivatives are with respect to $z$ and $\mathit{C}$ is the jet interfacial curvature, expressed as,

[TABLE]

For the fixed nozzle inlet, equation (5) is subject to boundary condition $h_{b}=1$ at $z=0$. Two more boundary conditions are needed to well define this differential problem of order three. However, exempting the jet tip from the base flow calculation gives us the liberty to impose a constant slope ($h_{b}^{\prime}=0$) and curvature ($h_{b}^{\prime\prime}=0$) at the exit of the jet. It should be noted that the boundary conditions applied at the jet exit should be treated as a way to close the differential problem rather than depicting physical boundary conditions. We made sure that these boundary conditions did not impact the overall base state solution by computing the solution over a large enough domain where the base state solution naturally converges to a solution with $h_{b}^{\prime}=0$ and $h_{b}^{\prime\prime}=0$.

2 Nonlinear simulations

The strength of a nonlinear simulation lies in its ability of capturing the exact response of the jet interface, including the shape close to the breakup point where the interface height $h$ approaches to a zero value. Often, the external forcing does not result in the breakup of fixed sized drops, rather the regular sized drops are followed by much smaller ‘satellite drops’.

Keeping this in view, we analyze the response of the jet in presence of an external forcing. We aim at finding the optimal forcing which results in the most unstable jet. The breakup length, which is the length of the stable jet between the nozzle and the breakup point, is chosen as the quantifier to compare the effect of different forcing, with the optimal forcing resulting in the shortest possible breakup length.

We begin with the description of the modified non-linear governing equations used for the simulations followed by the numerical scheme implemented to capture jet breakup. Finally we present the comparison of breakup characteristics of the jet for different inlet forcing .

2.1 Governing equations

In order to remove the singularity in expression (2) for the pressure, when $h(z,t)\to 0$, we define the interface height $h(z,t)$ in terms of function $a(z,t)$ where $a=h^{2}$. The governing equations (1) thus transform into,

[TABLE]

The base state solution for the jet interface $(a_{b}=h_{b}^{2})$ is obtained by solving equation (5). We model the external forcing on the jet by perturbing only the inlet velocity using a forcing of the form,

[TABLE]

where $\epsilon$ represents the amplitude of the forcing and $\omega$ represents the angular forcing frequency. In presence of the forcing, the boundary conditions at the inlet are modified to $a(0,t)=1$ and $u(0,t)=\sqrt{\mathit{We_{in}}}+u_{f}$. No boundary conditions are defined at the other extremity of the domain close to the tip. Nonetheless, a special treatment is applied for the tip as explained in the next section.

2.2 Numerical scheme

The governing equations (7) are first discretised in space, after which the resulting ordinary differential equations (ODE) are integrated in time. Diffusion terms are evaluated using second-order finite differences, with a central scheme for intermediate nodes and a forward or backward scheme for boundary nodes. Advection terms are obtained using a weighted upwind scheme inspired by Spalding (1972) hybrid difference scheme. Unlike the latter, which approximates the convective derivative using a combination of central and upwind schemes, we evaluate the derivative based on a combination of forward and backward finite differences. An advection term $dA/dz$ is evaluated at node $i$ as

[TABLE]

where indices $b$ and $f$ refer to the backward and forward finite difference schemes, and $\beta$ is a weight coefficient that depends on the local value of velocity $u$ at node $i$ together with a parameter $\alpha$,

[TABLE]

For the range of feed velocities considered in this study, numerical stability was always ensured by using a 10-point stencil. Thus, the backward difference term relies on a stencil that spans nodes $i-5$ to $i+4$, and the forward difference term employs nodes $i-4$ to $i+5$. For large enough downstream or upstream velocities, $\beta$ will tend to $1$ or [math] respectively; hence (9) reduces to a regular upwind difference scheme. For smaller velocity magnitudes in between, (9) produces a weighted combination of backward and forward differences. In our simulations, we choose $\alpha=50$ so that the transition between the backward and forward difference schemes mostly occurs when $|u|<0.05$. Finally, advection terms at nodes close to the boundary are evaluated based on the values of the closest 9 adjoining nodes.

After obtaining all spatial derivatives, the resulting ODEs are integrated using the MATLAB solver ode23tb, which implements a trapezoidal rule and backward differentiation formula known as TR-BDF2 (Bank et al., 1985), and uses a variable time step to reduce the overall simulation time.

The numerical domain $L$ is taken sufficiently large to capture the breakup of the jet. The jet interface is initialised by the solution of (5) obtained numerically with the MATLAB bvp4c solver. The validation of the numerically obtained base state solution is presented in Appendix 7.1. It should be noted that the steady state is implemented only for a part of the numerical domain and the interface is initialised to 0 for the remaining part. The axial span of the base state solution does not affect the quasi-steady jet characteristics, which are the focal point of our numerical analysis. A validation for the same is presented in Appendix 7.2.

At every time step, the solution is evaluated for three conditions: (i) Pinch-off (breakup): It is defined as when the value of $a$ passes below a threshold value of $10^{-5}$. The corresponding time $\mathit{t_{po}}$ is saved and the position of the jet tip is updated as $\mathit{N_{tip}=N_{po}}$, where $\mathit{N_{po}}$ is the pinch-off location. The solution for $a$ and $u$ beyond $\mathit{N_{tip}}$ is set to zero. For subsequent time steps, $\mathit{N_{tip}}$ has two possibilities – it can either advance or recede, which requires the following two conditions. (ii) Advancing jet: The values of $a$ at nodes $\mathit{N_{tip}-1}$ and $\mathit{N_{tip}}$ are extrapolated to find $a$ at $\mathit{N_{tip}+1}$. If the extrapolated value is larger than a predefined value of $5\cdot 10^{-3}$, the parameter $\mathit{N_{tip}}$ is incremented by 1, and $a$ and $u$ at the new $\mathit{N_{tip}}$ are assigned values extrapolated from its previous two neighbours. (iii) Receding jet: If the value of $f$ at $\mathit{N_{tip}}$ falls below a predefined value of $10^{-3}$, $a$ and $u$ at $\mathit{N_{tip}}$ are set to zero and the parameter $\mathit{N_{tip}}$ is reduced by 1. These three conditions enable the numerical integration of the governing equations in a way that captures accurately the breakup of the jet and the motion of the tip. A validation of the numerical scheme is presented in Appendix 7.3.

2.3 Nonlinear simulation results

Using the numerical scheme presented in the previous section, nonlinear simulations were performed for a jet governed by equation (7) for fixed inlet characteristics: $\mathit{Oh_{in}}=0.3$, $\mathit{We_{in}}=1.75$ and $\mathit{Bo_{in}}=0.1$. The jet inlet velocity is subjected to time harmonic forcing of the form given by equation (8) with a fixed amplitude $\epsilon$ and for forcing frequency $\omega=[0.4-3.2]$

The simulations were run for a sufficiently long time to enter a permanent regime wherein the jet breaks up at regular intervals of time and at fixed axial location. In the quasi-steady regime, the breakup period $\mathit{\Delta T_{po}}$ is defined as the time difference between two consecutive breakups or pinch-offs and the breakup length $l_{c}$ as the stable length of the jet between the nozzle and the pinch-off location. We use $l_{c}$ as the quantifier to determine the stability of the jet to external forcing such that the most amplified disturbance caused by the optimal frequency $\mathit{\omega_{opt}}$ will compel the jet to have the shortest possible breakup length.

We begin our analysis for fixed amplitudes $\epsilon=10^{-2}$ and $\epsilon=10^{-4}$. The response of the jet due to different forcing frequency $\omega$ in the permanent regime for the two above mentioned amplitudes can be seen in Fig. 1. Jets enforced by the same disturbance amplitude at the inlet are represented by the same colour. For visual clarity we plot the response only for certain frequencies and for the jet shape pertaining to the shortest breakup length. First, Fig. 1 clearly shows the existence of a main drop and a satellite drop for all the frequencies. Second, for $\epsilon=10^{-2}$ we conclude that the optimal forcing frequency is $\mathit{\mathit{\omega_{opt}}}=1.38$ because it manifests the jet to have the shortest $l_{c}$. Third, and most strikingly, we notice that for a lower forcing amplitude of $\epsilon=10^{-4}$, the optimal forcing increases to $\mathit{\omega_{opt}}=1.68$. Finally, at all forcing frequencies, the breakup length for jet with $\epsilon=10^{-4}$ is always larger than for $\epsilon=10^{-2}$.

To investigate further the breakup characteristics for the amplitudes $\epsilon=10^{-2}$ and $\epsilon=10^{-4}$ due to $\mathit{\omega_{opt}}$, we plot the interface evolution in the permanent regime as shown in Fig. 2(a) and Fig. 2(c), where regular sized main drop formation is followed by the release of a satellite drop. For both the amplitudes, we see a distinct difference between the main and satellite drop radius.

The breakup period $\mathit{\Delta T_{po}}$ resulting from the forcing imposed in Fig. 2a and 2c are plotted in Fig. 2b and 2d, respectively, where the black triangle represent the $\mathit{\Delta T_{po}}$ obtained for two consecutive pinch-offs whereas the blue circles denote $\mathit{\Delta T_{po}}$ obtained for two consecutive pinch-offs of the same group, that is to say between two consecutive main (or satellite) drops. From the figure we conclude that even though the breakup period is the same for the group of main and satellite drops (as shown by blue circles) the time of formation of a satellite drop does not lie exactly midway between the time of formation of the main drops and vice versa. This results in obtaining two oscillating breakup periods (as shown by the black triangles). We further observe that the frequency of breakup ($\omega_{po}=2\pi/\Delta T_{po}$) obtained using the breakup period for consecutive main (or satellite) drops responds to the externally applied forcing at the jet inlet with $\omega_{po}=1.38$ and $1.68$ for $\epsilon=10^{-2}$ and $10^{-4}$, respectively.

Finally, for the constant flow rate of the jet, the breakup period related to the consecutive pinch-off’s is used for obtaining the drop radius for the satellite and main drops. We notice that at the optimal forcing frequency, the main drop radius $R_{m}$ decreases from 1.62 to 1.48 dimensionless units as $\epsilon$ reduces from $10^{-2}$ to $10^{-4}$. On the contrary, the satellite drop radius $R_{s}$ increases from 0.78 to 1.08 dimensionless units for $\epsilon=10^{-2}$ and $\epsilon=10^{-4}$ respectively. The longer intact jet length obtained for lower forcing amplitude $\epsilon=10^{-4}$ results in a larger downstream velocity close to the tip due to the presence of gravity. Eventually, it results in the formations of highly stretched satellite drops in comparison to the ones obtained for lower amplitude of $\epsilon=10^{-2}$ as seen in Fig. 2a and 2c.

We now return to the most salient feature observed in Fig. 1, where the optimal forcing frequency $\mathit{\omega_{opt}}$ increased with a decrease in forcing amplitude. To explore if this effect existed for smaller amplitudes, we simulated the same system for different forcing amplitudes $\epsilon=[10^{-2}-10^{-8}]$, and plotted the breakup length $l_{c}$ as a function of the forcing frequency $\omega$ as shown in Fig. 3, where the optimal forcing frequency for a fixed $\epsilon$, is marked with a red cross. The results show an increase in $l_{c}$ and $\mathit{\omega_{opt}}$ as $\epsilon$ decreases. The increase in breakup length is obvious due to the decreasing destabilizing strength of the forcing amplitude. The increase in optimal forcing frequency, however, is the most interesting observation drawn from the numerical results, since it is expected to saturate for small enough forcing amplitudes. We believe that the increase in $\mathit{\omega_{opt}}$ as $\epsilon$ decreases from $10^{-2}$ to $10^{-8}$ is a consequence of the stretched base state due to gravity, which results in the downstream stretching of the perturbation wavelength initiated at the nozzle. As the forcing amplitude decreases, the stable jet length $l_{c}$ increases and so does the stretching close to the jet tip. Thus to compensate for the larger stretching, the breakup potential of the forcing is sustained by increasing the forcing frequency.

To conclude, the numerical simulations confirm the dependence of the $\mathit{\omega_{opt}}$ on the forcing amplitude, a factor generally neglected for linear stability analysis as long as $\epsilon\ll 1$. The trend also constitutes a major difference from a jet with no gravity effect ($\mathit{Bo_{in}}=0$), where $\mathit{\omega_{opt}}$ is independent of $\epsilon$ (refer Fig. 20 in Appendix 7.4. Finally, we confirm that the preferred-mode analysis carried out for the jet is solely due to the effect of external forcing. For the given parameter range, the jet tip did not induce any self-sustained breakups.

3 Local stability analysis

The linear stability theory is applicable for small forcing amplitudes ($\epsilon\ll 1$) and does not take into account its absolute value, a parameter that has already been shown in Section 2.3 to influence the optimal forcing frequency. Nevertheless, we begin our analysis in the local framework by first obtaining the dispersion relation for parallel viscous jets in absence of gravity, which is used as a basis for obtaining the absolute/convective instability transition criteria in Section 3.2. Next, the dispersion relation for parallel jets is suitably modified to include the spatial variation of the gravitationally stretched base flow and the spatial stability of the jet is performed in Section 3.3. Since the base state is spatially evolving, we extend our stability analysis using the WKBJ formulation in Section 3.4.

3.1 Local stability analysis for jets in absence of gravity

We derive the dispersion relation for the coupled equations (1), governing the growth of small perturbations about the base state. Considering the normal mode expansion, the flow variables $h(z,t)$ and $u(z,t)$ are decomposed as:

[TABLE]

where $\epsilon\ll 1$, with $\hat{h}$ and $\hat{u}$ as complex constants. $k$ and $\omega$ are respectively the dimensionless spatial wavenumber and the temporal frequency, which may both be complex. Similarly, for the variable representing the square of interface, $a(z,t)$ is decomposed as,

[TABLE]

where $a_{b}=h_{b}^{2}$ and $\hat{a}=2h_{b}\hat{h}$. Inserting the above expansion into equation (1), linearizing about $(h_{b},u_{b})$, and replacing $h^{2}\rightarrow a$, will lead to a linearised system of equations which can be formulated as an eigenvalue problem with the eigenmodes represented by $\mathbf{{\hat{q}}}=[\hat{a},\hat{u}]$.

In absence of gravity, $\mathbf{{\hat{q}}}$ represents a constant independent of $z$, $u_{b}=\sqrt{\mathit{We}}$ and $h_{b}=1$. Combining both the linearised continuity and momentum equations leads to the dispersion relation,

[TABLE]

where $\mathit{Oh}$ is constant throughout the domain. The dispersion relation is used to perform a spatio-temporal stability analysis which includes the effect of advection speed of the jet on its stability properties. In this framework, we define the impulse response of a system to a localised perturbation which generates a wave packet growing in space and time. In the laboratory framework, the spatio-temporal behaviour of the wave packet can be described in terms of the complex absolute wave number $k_{0}$ and the corresponding complex absolute frequency $\omega_{0}=\omega(k_{0})$ whose imaginary part $\omega_{0,i}$ will determine the temporal evolution of the wave packet. For $\omega_{0,i}>0$ the system is absolutely unstable since the disturbance grows fast enough to invade entire domain in the laboratory frame and for $\omega_{0,i}<0$ the system is convectively unstable as the localised perturbations are allowed to convect downstream before they grow in the laboratory framework. The complex pair ($k_{0},\omega_{0}$) is defined using the saddle point condition or the Briggs-Bers zero-group velocity criterion, together with the dispersion relation

[TABLE]

where $\Delta$ represents the dispersion relation (13) and

[TABLE]

Equation (14) identifies the critical dimensionless speed $\mathit{We_{crit}}$, for a fixed $\mathit{Oh}$, which signifies the transition of the jet from an absolutely unstable to convectively unstable system as shown in Fig. 4 with the full and dashed lines. The two curves are obtained for a system initialised either using a low or a high $\mathit{Oh}$, respectively. For the intermediate values of $\mathit{Oh}$ we obtain two saddle points thus giving two distinct values of the critical curve.

3.2 Local stability analysis for jets in presence of gravity

Extending the formalism for parallel jets to spatially varying jets, we derive the dispersion relation for the coupled equations (1), governing the growth of small perturbations about the base state. The linearised system of equations obtained around the spatially varying base flow, can be formulated as an eigenvalue problem with the eigenmodes represented by $\mathbf{{\hat{q}(z)}}=[\hat{a}(z),\hat{u}(z)]$ as functions of $z$. Next we express the local stability of a gravity jet by introducing the terms $\mathit{Oh_{z}}$ and $\mathit{\mathit{We_{z}}}$, which are the local dimensionless numbers at an axial distance of $z$ from the nozzle. They are expressed as,

[TABLE]

We then plot the values $\mathit{Oh_{z}}$ and $\mathit{We_{z}}$ along the entire axial domain $L$ above the absolute-convective transition curve in Fig. 4. The variation in local $\mathit{Oh_{z}}$ and $\mathit{\mathit{We_{z}}}$ along the jet defined within a domain size $L=50$, for $\mathit{Oh_{in}}=0.3$ and $\mathit{Bo_{in}}=0.1$ and for three different inlet Weber numbers, $\mathit{We_{in}}=[1.75,0.25,0.002]$ are represented by the markers in Fig. 4. The gap between consecutive markers is representative of an axial interval of 10 units. In each case, the red cross represents the inlet of the jet whose base state is shown in Fig. 5. We remind the reader that the case with $\mathit{We_{in}}=1.75$ corresponds to the jet whose numerical analysis has been presented in Section 2.3.

For $\mathit{We_{in}}=1.75$ and $0.002$, the entire jet exists in the convective and absolute region, respectively. For intermediate $\mathit{We_{in}}=0.25$, there exists a small pocket of absolute instability close to the nozzle, after which the local parameters modify along the downstream direction resulting in the transfer of the jet into a convectively unstable regime.

The parameter $\mathit{Bo_{in}}$ indirectly decides the instability of the jet by affecting the base state solution. Since $\mathit{Bo_{in}}$ is constant, its relative strength for the stretching of the jet interface depends on the corresponding value of $\mathit{We_{in}}$, with the effect being more pronounced for lower values of $\mathit{We_{in}}$ as shown in Fig. 5.

Next, for the spatially varying base flow, we perform the stability analysis in a local framework wherein the system is considered parallel at each axial location. The local dispersion relation, which now includes the local spatially varying base flow properties and is given by,

[TABLE]

For the convectively unstable jet ($\mathit{We_{in}}=1.75$), the solution of the dispersion relation (17) for a given range of complex $\omega$ (with $\omega_{i}>0$) results in obtaining four spatial branches which are expressed as the roots of the fourth-order polynomial (17). The solution consists of upstream (referred as $k^{-}$) and downstream (denoted by $k^{+}$) propagating branches. To identify these branches, we successively add an artificial $\omega_{i}$ so as to separate the branches into the upper $k_{i}>0$ and lower $k_{i}<0$ planes (Huerre & Rossi, 1998; Gallaire & Brun, 2017). For a downstream propagating $k^{+}$ branch damped in space, the associated $k_{i}>0$. Based on this analysis, we obtain two downstream and two upstream propagating waves for the dispersion relation (17). The $k$ branches for the localised dimensionless numbers at the nozzle inlet ($z=0$) and domain end ($z=50$) are shown in Fig. 6(a) and 6(b) respectively with the two $k^{+}$ branches denoted by the black and green colour and the two $k^{-}$ waves by the red and blue colour. The presence of two $k^{+}$ and $k^{-}$ waves is not specific to the present jet characteristics but rather exists for all the tested cases in the range of $\mathit{Oh_{in}}=[0.1~{}10]$, $\mathit{We_{in}}=[0.8~{}10]$, $\mathit{Bo_{in}}=[0~{}1]$ for $L=50$.

3.3 Spatial stability analysis

Since the base flow with $\mathit{We_{in}}=1.75$ exists in the convectively unstable regime (see Fig. 4), we then proceed to analyse the base flow using the spatial stability framework, wherein the spatial growth rate for the imposed real frequency determines the flow stability.

Given the polynomial nature of the dispersion relation, there are four spatial waves. We have verified that two of them are $k^{+}$, downstream propagating, waves while the remaining two are $k^{-}$, upstream propagating, waves (see Fig. 6). For a more detailed account on the nature of spatial waves in capillary jets, depending on the flow model, the reader is referred to Guerrero et al. (2016).

Among the four $k$ waves, only one of the $k^{+}$ waves is seen to be amplified. To obtain this dominant $k$ wave we plot the spatial growth rate $-k_{i}$ as a function of the real forcing frequency $\omega$ at nozzle inlet. As shown in Fig. 7(a), among the four $k$ branches, only the branch denoted in black has a growth rate which is positive in its propagation direction. We chose the $k$ wave corresponding to this amplified $k^{+}$ branch as the dominant wavenumber for all frequencies. The relevant $k(z)$ branches are then obtained for different $z$ along the jet as shown in Fig. 7(b) and 7(d) by imposing the spatially dependent base flow and $\mathit{Oh_{z}}$ in equation (17). Fig. 7(b) shows that the most amplified frequency shifts to higher values as one travels away from the nozzle. The associated eigenmode $\hat{\bf q}(z)$ also changes as one progresses downstream. Imposing $\left\lVert\mathbf{\hat{q}}\right\rVert=1$ at every axial location as the normalisation condition, together with $\hat{a}_{i}=0$ to set the phase, we see in Fig. 7(c) the evolution of the locus of the real and imaginary parts of the remaining degrees of freedom $\hat{u}_{r}$ and $\hat{u}_{i}$ as $z$ increases (remember that $\hat{u}_{r}^{2}+\hat{u}_{i}^{2}+\hat{a}_{r}^{2}=1$). While this locus is difficult to interpret from a physical point of view, it highlights the change of the eigenmode along the jet axis in such nonparallel gravity driven jets.

The knowledge of the relevant $k^{+}$ wave obtained for a given $\omega$ allows us to evaluate the leading order response due to different forcing frequencies imposed on the base flow, conveniently expressed as

[TABLE]

The overall response norm defined in a domain size $L$ is then given as,

[TABLE]

This allows us to determine the optimal forcing frequency $\mathit{\omega_{opt}}$ which results in the maximal gain,

[TABLE]

attained at a frequency $\omega_{opt}$. Fig. 8 (in dotted lines) shows the spatial gain as a function of forcing frequency for two arbitrary domain sizes $L=50$ and $60$. We notice that $\mathit{\omega_{opt}}$ shifts from 1.16 to 1.21 as we increase the domain size.

3.4 Weakly nonparallel stability analysis (WKBJ)

In order to further incorporate the non-parallelism of the base flow, we extend our spatial analysis by including the WKBJ formalism introduced by Gaster et al. (1985) and Huerre & Rossi (1998) for a spatial mixing layer and applied by Viola et al. (2016) for swirling flows.

In this framework, we introduce a slow streamwise scale $Z$, which relates to the fast scale $z$ as $Z=\eta z$, where $\eta\ll 1$ is a measure of the weak non-parallelism. The new base flow depends only on $Z$ and the global response to inlet forcing takes the modulated wave form:

[TABLE]

where $\mathbf{\hat{q}}(\omega,Z)$ is the local eigenmode and $k(\omega,Z)$ the local wavenumber at section $Z$ and a fixed forcing frequency $\omega$. The amplitude function $A(Z)$ acts as an envelope, smoothly connecting the progressive slices of the parallel spatial analysis. At each axial location, we impose $\mathbf{\hat{q}}^{H}\cdot\mathbf{\hat{q}}=1$, where $(\cdot)^{H}$ is the transconjugate. As described in Appendix 7.7, imposing an asymptotic expansion and a compatibility condition, the local stability analysis is retrieved at zeroth-order in $\eta$, while at first order in $\eta$ the following amplitude equation is obtained:

[TABLE]

whose solution is given as

[TABLE]

The functions $M(Z)$ and $N(Z)$ are defined in Appendix 7.7. The amplitude at the inlet is set as, $A(0)=1$ which simplifies the forcing expression at the inlet to $\mathbf{q}^{\prime}(0,t)=\mathbf{\hat{q}}(0)\exp(i\omega t)$. Finally we express the total spatial gain at first order as

[TABLE]

The global gain of the response due to the forcing frequency, for fixed domain sizes, is reported in Fig. 8, where $\mathit{\omega_{opt}}=1.24$ and $1.30$ for $L=50$ and $60$, respectively. The WKBJ approximation greatly modifies the gain when compared to the zeroth-order analysis and shifts the optimal forcing frequency predicted from the spatial analysis which excludes the amplitude equation. However, to truly assess the validity of the amplitude equation one needs to analyse the base flow in the global framework using the resolvent analysis, which will be the focus of our discussion in the next section.

4 Global stability analysis

Unlike the local stability analysis, the global stability framework allows taking into consideration the axially varying base state due to the stretching effect of gravity. In this framework, we first evaluate the inherent global stability of the base flow in Section 4.1. We next perform a resolvent analysis in Section 4.2 on the globally stable base flow to evaluate its response in presence of a given perturbation.

4.1 Global stability

Since the base flow is spatially varying, the perturbations imposed on it are no longer sought in the form of Fourier modes but are expanded in the form:

[TABLE]

where $\epsilon\ll 1$ and $\tilde{h}(z),\tilde{u}(z)$ are the global stability modes related to the complex growth-rate $\lambda$. Substituting expressions (25) in equations (1) and linearising around the base state $(h_{b},u_{b})$ results in the general eigenvalue problem of the form

[TABLE]

with boundary conditions $\tilde{h}(0,t)=0$ and $\tilde{u}(0,t)=0$. We do not impose any boundary conditions at the end of domain $z=L$ since it is not possible a priori to distinguish between amplifying perturbations and transient disturbances. This will occur in any problem that involves an ‘active system’ and can support amplifying waves (Briggs (1964) and Leib & Goldstein (1986)). The global stability analysis presented in Rubio-Rubio et al. (2013) does not impose any boundary conditions for $z=L$ since the numerical method naturally converges to the most regular asymptotic solution of the base flow equation (5) and eigenvalue problem (26), as $z\to\infty$. Nonetheless, we checked that the dominant eigenvalue and eigenmode were unaffected by the presence of Neumann boundary condition at $z=L$, namely $\frac{d\tilde{u}}{dz}(L)=0$.

The complete expressions for the linear operator $\mathit{M}$ can be found in Appendix 7.5. The solution of equation (26) results into a set of eigenmodes ($\tilde{h},\tilde{u}$), whose growth rate and frequency are given by the real ($\lambda_{r}$) and imaginary ($\lambda_{i}$) parts of the related eigenvalue. A base state is stable to self induced oscillations provided $\lambda_{r}<0$.

To solve the eigenvalue problem, the Chebyshev collocation method is used for obtaining the differential operators. Derivatives with respect to z are calculated using the standard Chebyshev differentiation matrices. Denoting the non dimensional physical domain as $L$, the domain is mapped into the interval $-1\leqslant y\leqslant 1$ by using the transformation $z=[(L/2)\times(y+1)]$. A validation of the global scheme with the results of Rubio-Rubio et al. (2013) is presented in Appendix 7.6.

For the three cases of jets described in Fig. 5, with $L=50$, $\mathit{Oh_{in}}=0.3$, $\mathit{Bo_{in}}=0.1$ and three different values of $\mathit{We_{in}}$ a global stability analysis is carried out using different resolutions ($N_{1}=100,N_{2}=125,N_{3}=150$) to exclude spurious eigenvalues. The eigenvalue spectrum for $\mathit{We_{in}}=0.25,1.75~{}\text{and}~{}0.002$, are represented in Fig. 9(a), 9(c) and 10 respectively. The dominant eigenvalues have $\lambda_{r}<0$ (for $\mathit{We_{in}}=0.25$ and $1.75$) and $\lambda_{r}>0$ (for $\mathit{We_{in}}=0.0025$) thus representing globally stable and unstable jets, respectively. Note however that the local stability analysis of the globally stable flow with $\mathit{We_{in}}=0.25$ predicts the jet to have a small ‘pocket’ of absolute instability close to the nozzle.

Eigenmodes corresponding to the dominant eigenvalues are presented in the accompanying figure. We note that the dominant eigenmode, as represented in 9(b), 9(d) and 10, has an amplitude that grows downstream. Fig. 10(b) also shows that the wavelength grows downstream. This is a consequence of the fluid acceleration caused by gravity (Tomotika, 1936 and Rubio-Rubio et al., 2013) and can be interpreted from Fig. 7(d) where $k_{r}$ is seen to decrease with increasing $z$ for a fixed forcing frequency $\omega$. Further we see that close to the outlet, the eigenmodes evolve at a much larger length scale compared to that related to the variations in steady state jet, thus strengthening the argument that weakly non parallel stability analysis should be used with care for predicting the global stability of the gravity jet.

4.2 Global resolvent

Analysing the linear response of the base state for an external harmonic forcing at frequency $\omega$ is only well defined if the linear operator is stable or in other words the base state is stable, where the imposed perturbations are allowed to travel downstream before spreading in the entire domain under consideration. Else the algebraically amplified solution is superimposed by the unforced naturally growing exponential mode. Keeping this in mind, in this section we present the resolvent analysis for the stable gravity jet ($\mathit{We_{in}}=1.75$). To compare our results with the nonlinear simulations of Section 2, we impose similar inlet forcing conditions and approximate the gain predicted by the resolvent analysis in terms of the forcing amplitude.

4.2.1 Problem formulation

The external force $f$ is modelled as an incoming perturbation in the form of an unsteady upstream boundary condition of the 1D Eggers & Dupont equation (1). The resulting linearised equation is represented as,

[TABLE]

where as in the eigenvalue problem (26) $\boldsymbol{s}=[h,u]$, $\rm I\!I$ represents the identity matrix, $\mathit{M}$ the linear operator (detailed in Appendix 7.5) and $\mathit{B_{f}}$ the operator which expresses the effect of the inlet forcing onto the bulk equation. Considering a time-harmonic forcing, $f=\tilde{f}\text{exp}(-\text{i}\omega t)$, results in an asymptotic flow response $\boldsymbol{s}=\tilde{\boldsymbol{s}}\text{exp}(-\text{i}\omega t)$ at the same frequency. Here $\tilde{\boldsymbol{s}}=[\tilde{h},\tilde{u}]$. Imposing these transformations in (27) we obtain,

[TABLE]

Equation (28) is subjected to two inlet and two outlet boundary conditions. The solution forced only in $u$ at the nozzle satisfies $\tilde{h}=0$ and $\tilde{u}=1$, while for the one that is forced in both $h$ and $u$, $\tilde{h}$ and $\tilde{u}$ can be chosen arbitrarily. We use the former when comparing the results with the nonlinear simulations (where the forcing was applied using the form (8)) whereas the latter when comparing to the spatial and WKBJ analysis of Section 3.3 and 3.4.

Based on the results of the local and global analysis at $\mathit{We_{in}}=1.75$, the convective instability of the flow ensures that at the outlet any existing $k^{+}$ branch, obtained from the spatial analysis, will be transmitted downstream. Since the relevant $k^{+}(\omega,L)$ branch for a given $\omega$ and at $z=L$ can be obtained from the local analysis of Section 3.3, we impose for the solution of the equation (28) at $z=L$ the spatial response, specifically, $\tilde{h}(L)=\hat{h}\exp(\text{i}kL)$ and $\tilde{u}(L)=\hat{u}\exp(\text{i}kL)$, $k$ being the unique root of the dispersion relation corresponding to a donwstream amplified wavenumber. It should be noted that for active systems, such as the jet falling under the influence of gravity, it is not possible a priori to impose unique boundary conditions at the outlet. Thus, substitution of the resolvent response by the spatial response at the outlet should be treated as an approach to close the differential problem of equation (28) rather than depicting the physical boundary conditions. To ensure that these boundary conditions do not affect the final response over a given domain size $L$, we impose them for a domain size $L^{\prime}>L$, such that the response for all the frequencies over $L$ is independent of the imposed boundary condition.

Finally we express the magnitude of the response $\tilde{\boldsymbol{s}}$ due to the externally applied forcing in terms of the gain $G$, with the maximum gain expressed as,

[TABLE]

attained at $\omega_{opt}$. To measure the amplitude of the response and the forcing, we define $\mathit{Q}$ and $\mathit{Q}_{f}$ as the weight matrices of the discretised energy norm ($\lVert{\tilde{\boldsymbol{s}}\rVert}^{2}=\boldsymbol{s}^{\dagger}\mathit{Q}\boldsymbol{s}$) and the forcing norm ($\lVert{\tilde{f}\rVert}^{2}=f^{\dagger}\mathit{Q}_{f}f$), respectively, obtained for the Chebyshev space on the physical domain $L$ which is mapped into the interval $-1\leqslant y\leqslant 1$ by using the transformation $z=[(L/2)\times(y+1)]$. $\mathit{Q_{f}}$ is a $2N\times 2N$ matrix enabling us to distinguish forcing on $u$ only or on both components $u$ and $h$. Following the optimization method using singular value decomposition (SVD) described in Marquet & Sipp (2010) and Garnaud et al. (2013), we then express the optimal gain using the following eigenvalue problem,

[TABLE]

whose leading eigenvalue solution $\lambda$ gives $G^{2}_{\mathit{max}}$ and the associated eigenmode solution $\tilde{f}$ yields the optimal normalised forcing amplitude in $(h,u)$ to be applied at the inlet.

Since the linear analysis is based on small perturbations, the exact amplitude of the perturbation is unaccounted for in the expression (29). For the resolvent analysis we then define $\mathit{G_{h,fu}}$ as the gain in ${h}$ from a solution forced only in $u$ and $\mathit{G_{q,fq}}$ as the gain in $\boldsymbol{{q}}$ for a solution forced in $\boldsymbol{{q}}$, where $\boldsymbol{{q}}=[{a},{u}]$. The two expressions for the gain: $\mathit{G_{h,fu}}$ and $\mathit{G_{q,fq}}$ are formulated to replicate the forcing and gain definitions in the nonlinear simulations (Section 2) and the spatial analysis (Section 3.3 and 3.4), respectively.

Looking for the gain in $\mathit{G_{q,fq}}$ requires the inclusion of additional operators $\mathit{P}$ and $\mathit{H}$ which express $\boldsymbol{{q}}$ in terms of $\boldsymbol{{s}}$, and are given by,

[TABLE]

The operator $\mathit{H}$ modifies the imposed boundary conditions in terms of $\tilde{a}$, whereas the operator $\mathit{P}$ adequately expresses the response in $\tilde{a}$ in terms of $\tilde{h}$ such that $\tilde{a}=2h_{b}\tilde{h}$. Additionally, to have an explicit comparison with the spatial analysis, we apply a forcing at inlet which is obtained as the eigenmode solution of the spatial problem in Section 3.3. Thus,

[TABLE]

The gain for the imposed forcing is then obtained as,

[TABLE]

Note however that even though this formalism allows a direct comparison with the spatial analysis, the gain $\mathit{G_{q,fq}}(\mathit{\omega})$ does not represent the maximum optimal gain since we do not impose the optimisation of the inlet forcing vector using an SVD formalism as was done in (30). Fig. 11 demonstrates the difference between the gain and $\mathit{\omega_{opt}}$ computed using the direct mode from spatial analysis (in black) and through an optimisation problem which solves for the optimal mode (in blue). Indeed the resolvent gain based on the optimised mode is much larger in magnitude.

4.2.2 Results: Comparison with spatial stability analysis

To replicate the type of forcing and the expression of gain used in the spatial analysis in Section 3.3 and 3.4, we impose the eigenmode solution at the nozzle exit obtained from the spatial analysis as the forcing vector in the resolvent analysis. The resulting gain $\mathit{G_{q,fq}}$ for two different fixed domain sizes $L=50$ and $L=60$ are shown in Fig. 8. We observe from the figure that the inclusion of the amplitude equation in evaluating the spatial response by far improves the estimation of the true gain obtained from the resolvent analysis. Moreover the predicted $\mathit{\omega_{opt}}$ producing the largest $\mathit{G_{q,fq}}$ from the WKBJ analysis is in close agreement with that of the resolvent analysis. The response norm obtained using the different approaches agrees qualitatively (see Fig. 22). Its non-monotonic behaviour at $\omega=1.5$ (see Fig. 22(c-d)) is well captured and is in accordance with the work presented in Lizzi (2016). However, the difference in gain between the three methods originates as a result of the quantitative disparity in the response obtained at different frequencies. The divergence between the spatial and resolvent analysis is due to the stretching effect of gravity on the base flow. For a parallel base flow, the results obtained from both the methods are found to be identical (as shown in Appendix 7.4, Fig. 21).

4.2.3 Results: Comparison with nonlinear simulations

We compare the resolvent analysis with the nonlinear simulations of Section 2.3. Classically, the optimal forcing frequency $\mathit{\omega_{opt}}$ resulting in the maximum gain, can be deduced by plotting the gain $\mathit{G_{h,fu}}$ as a function of $\omega$ for a fixed domain size $L$. For capillary jets however, the domain size over which the perturbation grows cannot be fixed a priori. It is merely an outcome of the analysis which should compare well with the value of $l_{c}$ measured in the nonlinear simulations.

In order to circumvent this lack of consistency and in absence of the knowledge of $l_{c}$, we first plot $\mathit{G_{h,fu}}$ as a function of increasing domain sizes $L$ and for fixed $\omega$ as shown in Fig. 12(a) where the gain $\mathit{G_{h,fu}}(\omega,L)$ is computed for $\omega=[1-2.5]$ with $\Delta\omega=0.01$ and for $L=[10-240]$ with $\Delta L=10$. This results in a bundle of constant frequency curves, intersecting each other at different locations in $L$. In Fig. 12(a) we now define the dominant frequency at a given $L$ as the frequency with the maximum gain at $L$. A close examination reveals that there is a continuous transition in the dominant frequency as one moves along increasing domain sizes. This is shown in Fig. 12(b) where for clarity we plot only the envelope $\mathit{G_{opt}}(L)$ of the dominant frequency for all values of $L$. $\mathit{G_{opt}}(L)$ is attained for $\mathit{\omega_{opt}}(L)$.

As discussed previously, $L$ represents the breakup location along the jet where the nonlinear effects appear. Broadly speaking, nonlinearity enters the system when a small perturbation $\epsilon$ gives rise to a response of the order of 1, which suggests to approximate $l_{c}$ by the value of $L$ at which

[TABLE]

In other words, the gain at the breakup location $L=l_{c}$ should be equal to $1/\epsilon$. Using equation (34), we locate the gain in Fig. 12(b) for different forcing amplitudes $\epsilon=[10^{-2}-10^{-6}]$. At the given value of $\mathit{G_{h,fu}}$, a horizontal projection on the dominant frequency envelope will then decide the optimal forcing frequency for the given $\epsilon$. Finally, a vertical projection on $L$ from the intersection point on the dominant frequency envelope will provide the relevant breakup length $l_{c}$ for the forcing amplitude $\epsilon$. Extracting the results from Fig. 12(b), we compare the optimal forcing frequency and the breakup length for different $\epsilon$ with the nonlinear solutions of Section 2.3 in Fig. 13. The close agreement between the two approaches shows the strength of the resolvent analysis in predicting the $\mathit{\omega_{opt}}$ and $l_{c}$ especially without any prior information from the nonlinear simulations. In Fig. 13 the small difference in values in the two methods can likely be attributed to ad-hoc definition of the required threshold for nonlinear effects to kick in and breakup to occur.

5 Response to white noise

Up to now, we were only interested in the response of the jet to an external disturbance characterised by a constant forcing frequency. However, in reality the external disturbance is more likely to be composed of a broadband frequency rather than being harmonic. Thus, to model this physical perturbation we carry out nonlinear simulations consistent with the scheme presented in Section 2 by exciting the jet at the nozzle by a white noise $\xi(t)$ defined in the time interval $[0~{}T]$ and formulated in a similar way as in Mantič-Lugo & Gallaire (2016). The white noise signal $\xi(t)$ is characterised by a constant power spectral density (PSD) $S_{\xi\xi}(\omega)=|\hat{\xi}(\omega)|^{2}$ where $\hat{\xi}(\omega)$ is the Fourier transform of $\xi(t)$ and has an infinite power $P$ defined as,

[TABLE]

where $\sigma$ is the variance. Even though a pure white noise has infinite power (as $S_{\xi\xi}(\omega)>0$), physical systems are usually characterised by a band-limited white noise. We thus filter the digital random signal $\xi_{d}(t)$ with a band limiting frequency $\omega_{b}/2\pi=1$ to obtain the band limited white noise $\xi_{b}(t)$ as shown in Fig. 14. For $\xi_{d}(t)$ the Nyquist frequency is set by $\omega_{N}/2\pi$ which depends on the time step ($\delta t$) of the signal, such that $\omega_{N}/2\pi=1/2\delta t$. Here we chose $\delta t=0.01$.

The noise $\xi_{b}(t)$ is normalised to have zero mean, unit variance and unit power, with a constant value for PSD, where $2|\hat{\xi_{b}}|^{2}=2\pi/\omega_{b}$ which completely depends on the band limiting frequency. Finally we impose this filtered white noise as an inlet velocity condition for the jet defined by $\mathit{Oh_{in}}=0.3$, $\mathit{Bo_{in}}=0.1$ and $\mathit{We_{in}}=1.75$ and governed by the equations (7) by replacing the boundary condition (8) with,

[TABLE]

where $\epsilon$ is the amplitude of the white noise signal. The forcing is applied at two different amplitudes $\epsilon=10^{-2}$ and $10^{-4}$ and for large times ($T=2000$) so as to achieve results which are time independent. For the MATLAB solver ode23tb with varying step size, the maximum time step size is set as $\delta t$ and white noise for intermediate time steps is obtained through interpolation.

At every pinch-off on the jet, we note the breakup length $l_{c}$, the pinch-off period $\mathit{\Delta T_{po}}$ and the drop radius $\mathit{R_{drop}}$ at the time of breakup. The distribution of the breakup characteristics is shown as as a histogram in Fig. 15 and 16 and compared with the expected response of the jet in presence of the pure $\mathit{\omega_{opt}(\epsilon)}$, which corresponds to $\mathit{\omega_{opt}}=1.38$ and $1.65$ for $\epsilon=10^{-2}$ and $10^{-4}$, respectively. The breakup characteristics for $\mathit{\omega_{opt}}$ have been discussed in Fig. 2 and are depicted by red bars in Fig. 15 and 16.

The drop size distribution shown in Fig. 15(a) highlights the two distribution peaks concentrated around $\approx 0.9$ and $\approx 1.65$, representing the group of satellite and main drops respectively. This behaviour also exists for smaller $\epsilon=10^{-4}$ where the radius is aggregated at $\approx 1.05$ and $\approx 1.45$. The results for the main drop size are coherent to the ones obtained in the presence of pure optimal forcing where $\mathit{R_{drop}}=1.62$ and $1.45$ for $\epsilon=10^{-2}$ and $10^{-4}$, respectively. Thus even in the presence of the white noise, the response of the jet is dominated by its expected behaviour at $\mathit{\omega_{opt}}$.

Unlike the drop radius, the peak of the distribution of breakup length obtained by imposing the white noise is not in close agreement with that of the optimal forcing as shown in Fig. 16(a) and (b). Yet, we clearly see that the distribution spectrum shifts to large values of breakup length as $\epsilon$ is decreased, a behaviour similar to the one predicted by $\mathit{\omega_{opt}}$ where $l_{c}$ increases from $\approx 50$ to $\approx 125$ as $\epsilon$ is decreased. Similar conclusions can be drawn for the comparison of $\mathit{\Delta T_{po}}$ between white noise forcing and forcing with $\mathit{\omega_{opt}}$ from 16(c) and (d) where we plot the breakup period between two consecutive drops.

6 Conclusion and perspectives

In this work, we inspect the response of a spatially varying gravitationally stretched jet subjected to an inlet velocity perturbation. The forcing is characterised through the frequency and the amplitude, the latter playing a major role in the determination of the optimal forcing frequency. The results of the numerical simulations performed on the nonlinear 1D Eggers & Dupont equations shows an increase in optimal forcing frequency and the breakup length as the forcing amplitude is decreased. We found that the amplitude dependent preferred mode is a characteristic of gravity driven jets only. A pure capillary jet, base state of which is independent of gravity-induced stretching, does not sustain such a behaviour. In such cases, decreasing the forcing amplitude only resulted in an increase of the breakup length with the optimal frequency remaining fixed at all amplitudes.

The linear stability theory characterised the jet flow used for nonlinear simulations as locally unstable and globally stable. Based on the absolute-convective transition criteria, we analysed the local stability at each section along the axial direction. The solution of the dispersion relation and the subsequent analysis for the downstream propagating spatial waves helped in confirming the predominant wave to be used for the zeroth order spatial gain expression. The strong non-parallelism of the base flow close to the nozzle motivated the incorporation of the WKBJ framework which markedly improved the prediction of the optimal forcing frequency in comparison to the resolvent analysis. However the spatial gain was still observed to be lower than the resolvent. As suggested by Le Dizès & Villermaux (2017), using advanced stability tools (Schmid, 2007) which accounts for non-parallel effects and non-modal growth leads to an estimation of a more realistic spatial response.

This task was tackled using a resolvent analysis which accurately captured the linear response of stable jets in presence of an external forcing. Assuming a simple global amplitude breakup threshold criterion, the linear resolvent analysis becomes capable in predicting both the breakup length and the optimal forcing frequency given the amplitude of the forcing. The results of the nonlinear simulations and the resolvent for different forcing amplitudes are quantitatively comparable, thus underlining the importance of the resolvent analysis. Besides forcing the jet inlet with a fixed frequency, we also studied the response to a white noise, to analyse its natural response to a distributed forcing frequency range. Surprisingly, even in the presence of the white noise, the dominant response of the jet is close to the one seen from the optimal frequency at that amplitude.

In presence of the external forcing, a dominant feature seen from the nonlinear simulations is the formation of a main and a satellite drop at the time of breakup. Nevertheless, to properly examine the consequence of the forcing amplitude on the final drop size, there is a need to enhance the nonlinear model by including the physics of drop coalescence and disintegration as done by Driessen & Jeurissen (2011). Post breakup, the state of the jet after the pinch-off should be inferred from the system before the breakup. Additionally, the choice for drop curvature is of paramount importance since a given breakup can possess variety of drops shapes-each on different length scales (Kowalewski, 1996).

On a different note, if the final aim is to eliminate the presence of satellite drops, the forcing should be modified such that it leads to the selective production of equisized drops. In this direction the work of Chaudhary & Redekopp (1980), who controlled satellite drops by forcing the jet with a suitable harmonic added to the fundamental; and Driessen et al. (2014), who controlled the size of the droplet breaking off from a parallel jet by imposing a superposition of two Rayleigh-Plateau-unstable modes on the jet, could serve as the basis for formulating a theory for the spatially varying gravity jets.

I.S. thanks the Swiss National Science Foundation (grant no. 200021-159957). The authors would like to thank Eunok Yim for extremely valuable discussions on the modelling of white noise disturbance and in the interpretation of the local/global response. The authors would also like to thank Tobias Ansaldi and Giorgio Rocca who worked on the foundation of the numerical code for parallel jets as well as Adrien Jean Pierre Bressy for efficiently performing several simulations on the enhanced version of the numerical code provided to him.

7 Appendix

7.1 Numerical base state solution validation

In this appendix, we show the validation of our numerically obtained base state solution of the governing equations (5) with the experimental results of Rubio-Rubio et al. (2013) for three different jet flows. The MATLAB bvp4c solver along with the boundary conditions stated in Section 1 accurately captures the stretching (necking) close to the nozzle due to the effect of $\mathit{Bo_{in}}$.

7.2 Effect of initial condition on breakup characteristics

This section demonstrates the effect on breakup characteristics due to different initial conditions of the jet. Using the scheme described in Section 2.2 we perform numerical solutions for a jet with $\mathit{Oh_{in}}=0.3$, $\mathit{We_{in}}=1.75$, and $\mathit{Bo_{in}}=0.1$ excited with a forcing of amplitude $\epsilon=10^{-2}$ and frequency $\omega=0.8$. In the first case, the jet is initialised as a circular tip of radius 1 (Fig. 18(a)) and in the second case with the base state solution obtained by solving equation (5) defined for an axial length of 100 (Fig. 18(b)). In both the cases the numerical domain is considered large enough to capture all the breakups. As shown in Fig. 18, both the jets with different initial conditions have different transient dynamics upto $t=55$ (tip) and $t=40$ (base state) after which they enter the permanent regime. In this regime, the breakup length $l_{c}$ and period $\mathit{\Delta T_{po}}$ are identical as shown in subplots (c) and (d), respectively. It is thus safe to conclude that in the permanent regime the jet breakup is independent of the initial base state solution.

7.3 Numerical scheme validation

In this appendix, we show the validation of our numerical scheme described in Section 2.2 for the simulations of reduced 1D Eggers & Dupont (1994) equations represented by the equation (7). For the purpose of validation, we use the numerical data of van Hoeve et al. (2010) which are described for micro-jets of initial radius $h_{0}=18.5\,\mu$m with density $\rho=1098\,\mathrm{kg}/\mathrm{m}^{3}$, viscosity $\eta=3.65\,$mPa.s, and surface tension $\gamma=67.9\,$mN/m. The jet is injected at a constant flow rate $Q=0.35\,$mL/min, corresponding to an initial jet velocity $U_{0}=Q/(\pi h_{0}^{2})=5.4\,$m/s. The flow can thus be described by the dimensionless numbers $\mathit{Oh_{in}}=0.1$ and $\mathit{We_{in}}=8.7$. Since the gravity effects are not considered in the experiment, we inject $\mathit{Bo_{in}}=0$ in equation (7).

To initiate jet breakup in their numerical simulations, a harmonic modulation of the dimensional nozzle radius is applied as follows:

[TABLE]

with $\delta/h_{0}\approx 0.005$ the forcing amplitude, and $n$ the driving frequency. The latter is selected to match the optimum wavelength $\mathit{\lambda_{opt}}$ for jet breakup, that is, $n=U_{0}/\mathit{\lambda_{opt}}$. To ensure a constant flow rate $Q$ through the nozzle, the dimensional velocity is modulated correspondingly as

[TABLE]

The amplitude of the wave imparted by the forcing at the nozzle grows until it equals the radius of the jet. Pinch-off or jet breakup is then defined as when the minimum width of the jet is below a predefined value set to $10^{-3}h_{0}$.

In our numerical simulations, we compute solutions to the governing equations (7) with the same harmonic forcing and flow parameters as in van Hoeve et al. (2010). A hemispherical droplet described by $h=({h_{0}}^{2}-z^{2})^{1/2}$ is used as initial condition for the shape of the jet, the tip of which is therefore initially at $z=h_{0}$. The velocity is initialised to $u_{0}$ everywhere along the jet. A fixed number of grid points, corresponding to a discretization size $dz=0.05$, is uniformly distributed throughout the entire domain. The final validation is presented in Fig. 19, which shows comparison of the time series of the dynamics of jet breakup obtained from our numerical scheme and the numerical results from van Hoeve et al. (2010).

For both figures, the evolution of the jet shape is shown at time intervals of $2\,\mu$s. Our numerical model predicts a breakup period of $25\,\mu$s and a breakup length of $856\,\mu$m. The results of van Hoeve et al. (2010), have a breakup period of about $26$ to $30\,\mu$s and a breakup length of about $800\,\mu$m. The error in breakup length between the two codes can be explained by the difference in grid size. Overall, Fig. 19 shows a good agreement between both results and validates our numerical scheme.

7.4 Comparison between resolvent and spatial analyses

In this section we briefly show the preferred forcing frequency of the jet discussed in Section 2.3 in the absence of gravity. The jet is characterised by $\mathit{Oh_{in}}=0.3$, $\mathit{We_{in}}=1.75$ and $\mathit{Bo_{in}}=0$. As shown in Fig. 20 nonlinear simulations for the governing equations (7) for the zero gravity case using different forcing amplitudes $\epsilon$ show that $\mathit{\omega_{opt}}$ is independent of the chosen forcing amplitude, a behaviour in contrast to the situation where gravity is present (previously shown in Fig. 3).

Fig. 21 shows the comparison of gain, in absence of gravity, as function of forcing frequency using spatial and resolvent analysis for two different domain sizes. For convenience we also plot the resolvent gain $\mathit{G_{h,fu}}$ expressed in terms of the forcing applied for the nonlinear simulations. We note that all the curves, irrespective of the domain size, predict the same optimal forcing frequency $\mathit{\omega_{opt}}=0.74-0.76$, a value close to the nonlinear prediction of Fig. 20. Moreover, unlike the situation with $Bo=0.1$, we notice that in the absence of gravity the magnitude of the gain at all frequencies is well captured by the spatial analysis.

7.5 Linear operator for eigenvalue problem

For the eigenvalue problem related to the global stability in Section 4.1, the matrix $\mathit{M}$ is expressed as

[TABLE]

where the expressions $\mathit{M}_{11}$, $\mathit{M}_{12}$, $\mathit{M}_{21}$ and $\mathit{M}_{22}$ denote the following differential equations:

[TABLE]

In the group of equations (40), I​I is the identity operator, $D^{n}\equiv d^{n}/dz^{n}$, $s(z)=[1+(h^{\prime}_{b})^{2}]^{-1/2}$ and

[TABLE]

7.6 Global stability validation

In this section, we present the validation of the numerical scheme used for the global stability analysis presented in Section 4.1. The validation is done against the results of Rubio-Rubio et al. (2013) where the stability analysis is based on the same 1D Eggers & Dupont (1994) equations but made dimensionless using different characteristic length and time scales. The results of Rubio-Rubio et al. (2013) are based on dimensionless numbers $\mathit{We_{in}}$, $\mathit{Bo_{in}}$ and Kapitza$(\Gamma)$. For the purpose of comparison we obtain the equivalent $\mathit{Oh_{in}}$ expressed as

[TABLE]

For the eigenvalue problem, a non dimensional domain length $L=120$ is considered. For obtaining the dominant eigenvalue, the solution was computed for different values of $N$ lying between $[150\,\,250]$. Figure 23 shows the validation of the eigenvalue spectrum with $N=200$, steady state and dominant eigenfunction for $\Gamma=5.83$ and $\mathit{Bo_{in}}=1.8$ and two different values of $\mathit{We_{in}}$. The results obtained from the present model are in good coherence with that obtained from Rubio-Rubio et al. (2013). Comparing Fig. 23(a) and (b), we can observe there is critical Weber number $\mathit{{We_{in}}_{c}}$, for which the jet becomes marginally unstable as the real part of the leading eigenvalue is slightly positive. For the given values of $\mathit{Bo_{in}}$ and $\mathit{Oh_{in}}$, the $\mathit{{We_{in}}_{c}}$ is therefore equal to $3\times 10^{-3}$. If the Ohnesorge number is fixed and the bond number is varied, we can get the corresponding $\mathit{{We_{in}}_{c}}$ for each value of $\mathit{Bo_{in}}$. Figure 24 represents the curve for the $\mathit{{We_{in}}_{c}}$, below which the jet becomes linearly unstable. The results are compared to those obtained from Rubio-Rubio et al. (2013).

7.7 WKBJ formulation for axisymmetric 1D Eggers & Dupont equations

7.7.1 Linearised equations

Considering linear perturbations ($a^{\prime},u^{\prime}$) in jet interface and velocity around the base flow ($a_{b},u_{b}$), the linearised system of equations is written as:

[TABLE]

where $D_{i}$ with $i=1..3$ are the differential operators with respect to $z$. Equation (43) can be reformed as

[TABLE]

where

[TABLE]

where the coefficients $A_{1..4}$ and $B_{1..7}$ are given as,

[TABLE]

In the above equation $S$ replaces the term $\frac{1}{4}{a_{b}}^{\prime}(z)^{2}+{a_{b}}(z)$.

7.7.2 Linearised equation expressed in terms of slow variable

For the WKB ananlysis, we then introduce the spatial scales. The fast spatial scale $z$ is replaced by the slow scale $Z$, such that $Z=\eta z$. The base flow is now expressed as a function of $Z$ such that $a_{b}(Z)$ and $u_{b}(Z)$. Let us consider the following normal mode expansion for the perturbation:

[TABLE]

Injecting the transformations (48a-d) into (43), the linearised equations on a weakly non-parallel baseflow equations are expressed through equations (49)-(50),

[TABLE]

where the continuity equation converts to:

[TABLE]

and the momentum equation transforms as,

[TABLE]

Defining $\mathbf{\hat{q}}^{(1)}=[\hat{a}^{(1)}~{}\hat{u}^{(1)}]$ and $\mathbf{\hat{q}}^{(2)}=[\hat{a}^{(2)}~{}\hat{u}^{(2)}]$, we now consider the asymptotic expansion:

[TABLE]

and inject it into the governing equations (49)-(50) to obtain the local stability problem at $\eta^{0}$ and $\eta^{1}$.

Order $\eta^{0}$ At zeroth-order in $\eta$, the local stability problem is retrieved:

[TABLE]

The system of equations represented by (52) can be reframed using the linear operator $\mathit{L}$, such that

[TABLE]

Substituting expression for $\hat{u}^{(1)}$ from (52)a into (52)b, we finally obtain:

[TABLE]

the solution of which gives the four roots of $k$ for a given $\omega$. The relevant $k$ branch is tracked as discussed in Section 3. For a given $\omega$ and a predetermined $k$, the solution of the linear problem (53) gives the response $\mathbf{\hat{q}}^{(1)}$, a parameter needed to solve the local stability problem at $\eta^{1}$.

Order $\eta^{1}$: At first order we obtain,

[TABLE]

where operator $\mathit{Q}$ can be split into two parts:

[TABLE]

Operator $\mathit{R}$ is expressed as:

[TABLE]

and $\mathit{S}$ is defined as,

[TABLE]

where the individual parameters are expressed as:

[TABLE]

As explained in Huerre & Rossi, 1998; Viola et al., 2016, in order to have solutions of the inhomogeneous equation $\mathit{L}[\mathbf{\hat{q}}^{(2)}]=\mathit{Q}[A\mathbf{\hat{q}}^{(1)}]$, the forcing term $Q$ should be in the image of the operator $L$. This implies that $Q$ should be orthogonal to the corresponding adjoint eigenfunction $\mathbf{\tilde{q}}^{(1)}$ of the adjoint operator $\tilde{L}$, with respect to the defined inner product,

[TABLE]

This leads to the amplitude equation,

[TABLE]

solving which we obtain the amplitude solution $A(Z)$ which should then be expressed in terms of the fast length scale $z$. Finally, at first order, the response if given by,

[TABLE]