2. Experimental set-up

The computer-aided design (CAD) model and the schematic of the experimental set-up are shown in figure 1. A fluid reservoir with a nearly constant pressure head delivers the fluid to the rotating nozzle via an intermediate chamber. A rotary swivel joint is used to connect the stationary fluid chamber to the rotating nozzle. The swivel joint is spun with a stepper motor through a timing belt at the desired speed. A circular plate connected to the piezo actuator is embedded in the fluid chamber for a well-controlled periodic perturbation of the nozzle exit velocity at a desired frequency. The jet is visualised with the help of a high-speed camera that is fixed in the lab frame of reference. The camera records at a frequency of 2500 frames per second to capture consecutive breakup events. The pixel size of the sensor is $11 \times 11 \,\mathrm {\mu }{\rm m}^2$ and the scale factor was found to be of $53\,\mathrm {\mu } {\rm m}\,{\rm px}^{-1}$, which corresponds to a magnification of 0.207. The exposure time of the camera is set to $75\,\mathrm {\mu }{\rm s}$, which limits the blurring of the jet to around 1.6 px. Since the Rossby number is close to unity, the blurring due to streamwise motion and the rotation of the jet is similar. Completely eliminating blurring was not possible so a compromise was made between blurring and good contrast. The jet is illuminated from the back by a LED panel. The velocity of the jet is determined at an estimated accuracy of $2\,\%$ from the drop in liquid level of the fluid reservoir. Properties of the working fluid are given in table 1.

Figure 1.

Experimental set-up. (a) Schematic of the set-up. (b) A CAD drawing of the piezo actuator and plate, fluid discharge tank and the rotary elbow. (c) Photo of the parts given in (b).

Table 1.

Experimental parameters and corresponding non-dimensional numbers used in the simulations.

3. Computational model

We performed a numerical study of the same flow geometry as considered in the experiments. The flow within the jet is decomposed into two components, namely, a stationary base state and perturbation. The perturbations are imposed at the nozzle and they evolve downstream along the jet. The jet trajectory and base flow are coupled with each other and, hence, must be solved simultaneously. Using the knowledge of the trajectory and how the jet velocity evolves downstream, we can then consider the perturbations on the jet. In what follows, we first present the mathematical framework in which we derive the governing equations in a moving frame of reference using curvilinear coordinates. Assuming a slender jet, we then derive the equations for the base flow and compute from this the fictitious forces acting on the flow along the jet trajectory. Finally, we derive the slender-jet model for the flow perturbations.

3.1. Mathematical framework

In a frame of reference co-rotating with the nozzle at an angular velocity $\boldsymbol {\varOmega } = (0,0,-\varOmega )$, the non-dimensional Navier–Stokes equations for the Newtonian liquid flow in the jet read

(3.1a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}$$

(3.1b)$$\begin{gather}\frac{{\rm D}\boldsymbol{u}}{{\rm D}t} ={-}\boldsymbol{\nabla} p + \frac{1}{Re}\boldsymbol{\nabla}\boldsymbol{\cdot}\left(2\boldsymbol{\mathsf{D}}\right) - \frac{1}{{\textit{Rb}}^2}\boldsymbol{\varOmega}\times\left(\boldsymbol{\varOmega}\times\boldsymbol{d}\right) - \frac{2}{{\textit{Rb}}} \boldsymbol{\varOmega}\times\boldsymbol{u} + \frac{1}{{\textit{Fr}}^2} \boldsymbol{g}, \end{gather}$$

where $\boldsymbol {g} = (0,0,-g)$ and $\boldsymbol{\mathsf{D}}=(\boldsymbol {\nabla }\boldsymbol {u}+\boldsymbol {\nabla }\boldsymbol {u}^{\rm T})/2$ is the rate-of-strain tensor and $\boldsymbol {d}$ is the distance vector with respect to the rotating axis. The definitions and values of the non-dimensional numbers are provided in table 1. The mathematical framework we use is very similar to what is described in detail in Shikhmurzaev & Sisoev (Reference Shikhmurzaev and Sisoev2017) and Noroozi et al. (Reference Noroozi, Arne, Larson and Taghavi2020), but briefly revisited here for completeness. Figure 2 shows a schematic of a curved jet. We start by describing a baseline (Shikhmurzaev & Sisoev Reference Shikhmurzaev and Sisoev2017), which is a spline that goes along the jet and remains inside it, using Cartesian basis vectors on a frame rotating with $\varOmega$ (see figure 2a),

(3.2)\begin{equation} \boldsymbol{D}(s) = X(s) \boldsymbol{\hat{x}} + Y(s) \boldsymbol{\hat{y}} + Z(s) \boldsymbol{\hat{z}}, \end{equation}

where $s$ is the arclength. From (3.2) the Frenet basis can be defined, i.e. a local orthogonal basis defined at every point on the baseline,

(3.3a)$$\begin{gather} \boldsymbol{T}(s) = \frac{{\rm d}\boldsymbol{D}}{{\rm d}s} = X'(s)\boldsymbol{\hat{x}}+Y'(s)\boldsymbol{\hat{y}}+Z'(s)\boldsymbol{\hat{z}} \quad \left(X'^2 +Y'^2 + Z'^2 = 1\right), \end{gather}$$

(3.3b)$$\begin{gather}\boldsymbol{N}(s) = \frac{{\rm d}\boldsymbol{T}}{{\rm d}s}\left\lvert\frac{{\rm d}\boldsymbol{T}}{{\rm d}s}\right\rvert^{{-}1}=\frac{X''\boldsymbol{\hat{x}}+Y''\boldsymbol{\hat{y}}+Z''\boldsymbol{\hat{z}}}{\sqrt{X''^2+Y''^2+Z''^2}}, \end{gather}$$

(3.3c)$$\begin{gather}\boldsymbol{B}(s) = \boldsymbol{T}\times\boldsymbol{N}, \end{gather}$$

where a prime denotes derivative with respect to $s$. The derivatives of the Frenet base vectors are related to each other by the Frenet formulas

(3.4)\begin{equation} \frac{{\rm d}}{{\rm d}s}\left[{\begin{array}{@{}c@{}} \boldsymbol{T} \\ \boldsymbol{N} \\ \boldsymbol{B} \end{array}}\right] = \left[{\begin{array}{@{}ccc@{}} 0 & \kappa & 0 \\ -\kappa & 0 & \tau \\ 0 & -\tau & 0 \end{array}}\right] \left[{\begin{array}{@{}c@{}} \boldsymbol{T} \\ \boldsymbol{N} \\ \boldsymbol{B} \end{array}}\right]. \end{equation}

Here $\kappa (s)$ and $\tau (s)$ are, respectively, the local curvature and the local torsion of the baseline, given by

(3.5a)$$\begin{gather} \kappa(s) = \sqrt{X''^2+Y''^2+Z''^2}, \end{gather}$$

(3.5b)$$\begin{gather}\tau(s) = \frac{X'(Y''Z'''-Z''Y''')+Y'(Z''X'''-X''Z''')+Z'(X''Y'''-Y''X''')}{(Y'Z''-Z'Y'')^2+(Z'X''-X'Z'')^2+(X'Y''-Y'X'')^2}. \end{gather}$$

Let us now define the curvilinear coordinate variables, $(s,r,\phi )$, where $s$ is again the arc length variable, and $(r,\phi )$ are the polar radius and angle, respectively, on the plane that is normal to the baseline at $\boldsymbol {D}(s)$. Note from the definition of the Frenet basis (3.3) that this is also the plane where $\boldsymbol {N}(s)$ and $\boldsymbol {B}(s)$ lie. So we can define an arbitrary point in the jet with a position vector measured from the origin in the Cartesian and curvilinear coordinates as follows:

(3.6a)$$\begin{gather} \boldsymbol{d}(x,y,z) = x\boldsymbol{\hat{x}} + y\boldsymbol{\hat{y}} +z\boldsymbol{\hat{z}}, \end{gather}$$

(3.6b)$$\begin{gather}\boldsymbol{d}(s,r,\phi) = \boldsymbol{D}(s) + r\cos{\phi}\boldsymbol{N}(s) +r\sin{\phi}\boldsymbol{B}(s). \end{gather}$$

Figure 2.

(a) Schematic of the top view of the set-up. (b) Close-up view of a section in the jet showing the different coordinate frames.

The free surface can be parametrically described as $r = R(s,\phi )$, with which the position vector to the free surface takes the form

(3.7)\begin{equation} \boldsymbol{d}(s,R(s,\phi),\phi) = \boldsymbol{D}(s) + R(s,\phi)\cos{\phi}\boldsymbol{N}(s) +R(s,\phi)\sin{\phi}\boldsymbol{B}(s). \end{equation}

On the free surface one has the balance of stresses in the normal and tangential directions,

(3.8a)$$\begin{gather} \boldsymbol{n}\boldsymbol{\cdot}\left({-}p+\frac{1}{Re}2\boldsymbol{\mathsf{D}}\right)\boldsymbol{\cdot} \boldsymbol{n} ={-}\frac{1}{{\textit{We}}}\kappa_{\textit{fs}}, \end{gather}$$

(3.8b)$$\begin{gather}\boldsymbol{t_i}\boldsymbol{\cdot}\left({-}p+\frac{1}{Re}2\boldsymbol{\mathsf{D}}\right)\boldsymbol{\cdot}\boldsymbol{n} = 0, \quad i = s,\phi, \end{gather}$$

where $\kappa _{\textit {fs}}$ is the curvature of the free surface. Here, the ambient phase is assumed to be dynamically inert and the effect of air drag is neglected.

We will use the notation $(s^1, s^2, s^3) = (s,r,\phi )$ to define the base vectors of the curvilinear coordinate frame as

(3.9)\begin{equation} \boldsymbol{e}_i = \frac{\partial \boldsymbol{d}}{\partial s^i} \quad (i=1,2,3). \end{equation}

Using (3.9), we can now express curvilinear base vectors in terms of the Frenet frame as follows:

(3.10)\begin{equation} \left[{\begin{array}{@{}c@{}} \boldsymbol{e}_1 \\ \boldsymbol{e}_2 \\ \boldsymbol{e}_3 \end{array}}\right] = \left[{\begin{array}{@{}ccc@{}} 1-r\kappa\cos{\phi} & -r\tau\sin{\phi} & r\tau\cos{\phi} \\ 0 & \cos{\phi} & \sin{\phi} \\ 0 & -r\sin{\phi} & r\cos{\phi} \end{array}}\right] \left[{\begin{array}{@{}c@{}} \boldsymbol{T} \\ \boldsymbol{N} \\ \boldsymbol{B} \end{array}}\right]. \end{equation}

The metric of this basis, $\boldsymbol{\mathsf{g}}$, defined as $g_{ij} = \boldsymbol {e_i \boldsymbol {\cdot} e_j}$, is given by

(3.11)\begin{equation} \boldsymbol{\mathsf{g}} = \left[{\begin{array}{@{}ccc@{}} (1-r\kappa\cos{\phi})^2 + (r\tau)^2 & 0 & r^2\tau \\ 0 & 1 & 0 \\ r^2\tau & 0 & r^2 \end{array}}\right]. \end{equation}

Note that the coordinate frame defined by (3.9) is not orthogonal as $g_{13}$ and $g_{31}$ are not zero in the presence of torsion or outside the baseline ($r\neq 0$). Another aspect of the curvilinear base is that two base vectors, namely $\boldsymbol {e}_1$ and $\boldsymbol {e}_3$, are not unit vectors since $g_{11}\neq 1$ and $g_{33}\neq 1$. As a consequence, there is a change in magnitude when a vector $\boldsymbol {v}=(v^1,v^2,v^3)$ is projected onto the curvilinear base. For instance, the projection of the velocity vector $\boldsymbol {u}$ on the curvilinear base is represented by $v^1\boldsymbol {e}_1+v^2\boldsymbol {e}_2+v^3\boldsymbol {e}_3 = u_s\boldsymbol {e}_1/|\boldsymbol {e}_1|+u_r\boldsymbol {e}_2/|\boldsymbol {e}_2|+u_{\phi }\boldsymbol {e}_3/|\boldsymbol {e}_3|$, where $u_s$, $u_r$ and $u_{\phi }$ are the physical components of the velocity using the normalized base vectors.

The inverse metric tensor ($g^{ij}$) is defined by $g^{ik}g_{kj}=\delta _j^i$, where $\delta _j^i$ is Kronecker delta, and is given by

(3.12)\begin{equation} \boldsymbol{\mathsf{g}}^{{-}1} = \frac{1}{\varDelta}\left[{\begin{array}{@{}ccc@{}} r^2 & 0 & -r^2\tau \\ 0 & \varDelta & 0 \\ -r^2\tau & 0 & (1-r\kappa\cos{\phi})^2 + (r\tau)^2 \end{array}}\right], \end{equation}

with $\varDelta = \det (g_{ij}) = r^{2} (r \kappa {(s )} \cos {(\phi )} - 1)^{2}$ the determinant of the metric tensor.

Using the metric and its inverse, we can express the gradient operator as (Wesseling Reference Wesseling2009)

(3.13)\begin{equation} \boldsymbol{\nabla} = g^{ij}\boldsymbol{e}_j\frac{\partial}{\partial s^i} = \left(\frac{1}{e_{1T}^2}\frac{\partial}{\partial s}-\frac{\tau}{e_{1T}^2}\frac{\partial}{\partial \phi}\right)\boldsymbol{e}_1 + \frac{\partial}{\partial r}\boldsymbol{e}_2 + \left(\frac{g_{11}}{r^2 e_{1T}^2}-\frac{\tau}{e_{1T}^2}\right)\boldsymbol{e}_3, \end{equation}

with $g_{11}=(1-r\kappa \cos {\phi })^2 + (r\tau )^2$ being the first element of the metric tensor given in (3.11) and $e_{1T}=1-r\kappa \cos {\phi }$ is the tangential component of the local base vector in the $s$ direction, namely the first element of the transformation matrix given in (3.10).

Finally, the coefficients to connect the variations of the base state vectors (also known as the Christoffel symbols of the second kind $\varGamma ^k_{ij}$, such that ${\partial \boldsymbol {e}_i}/{\partial s^j} = \varGamma ^k_{ij}\boldsymbol {e}_k$) can be written as

(3.14)\begin{equation} \varGamma^k_{ij} = \frac{g^{kl}}{2}\left(\frac{\partial g_{lj}}{\partial s^i}+\frac{\partial g_{li}}{\partial s^j}-\frac{\partial g_{ij}}{\partial s^l}\right). \end{equation}

The full expressions of each Christoffel symbol are provided in Appendix C. From here onwards, it is tedious book keeping to determine the terms in (3.1a) and (3.1b) in curvilinear coordinates. This was achieved by the help of the symbolic mathematics module SymPy (Meurer et al. Reference Meurer2017). Some of the terms in the conservation equations in general curvilinear coordinates are provided in Appendix C.

To pave the way for the slender-jet approximation, one needs to acknowledge the existence of two disparate length scales existing in the problem. These are namely the rotating arm, $L$, and the orifice radius, $R_0$. We scale $s$ with $L$ and $r$ with $R_0$. The velocity at the orifice, $U_0$, is used for scaling the velocity in the $s$ direction, so the advective time scale becomes $t_{\textit {adv}} = L/U_0$. Finally, the pressure is scaled with $\rho U_0^2$.

3.2. Jet trajectory and base flow

To obtain the leading order steady state equations for the spiralling jet trajectory, one needs to consider the limit where the jet slenderness $\epsilon$ tends to zero, i.e. $R_0/L\rightarrow 0$. The steady state jet base flow equations, namely, the mass conservation, the projections of the momentum conservation in tangential, normal and binormal directions and the arc length condition have the following form:

(3.15)\begin{gather} \frac{{\rm d}}{{\rm d}s}\left(UR^2\right) = 0, \end{gather}

(3.16)\begin{gather} U\frac{{\rm d}U}{{\rm d}s} ={-}\frac{1}{{\textit{We}}}\frac{{\rm d}}{{\rm d}s}\left(\frac{1}{R}\right) + \frac{3}{R^2Re}\frac{{\rm d}}{{\rm d}s}\left(R^2 \frac{{\rm d}U}{{\rm d}s}\right) \underbrace{-\frac{Z'}{{\textit{Fr}}^2}}_{f_{\textit{gravity}}} + \underbrace{\frac{XX'+YY'}{{\textit{Rb}}^2}}_{f_{\textit{centrifugal}}}, \end{gather}

(3.17)$$\begin{gather} \left(U^2-\frac{1}{RWe}-\frac{3}{Re}\frac{{\rm d}U}{{\rm d}s}\right)\left(X''^2+Y''^2+Z''^2\right) = \frac{XX''+YY''}{{\textit{Rb}}^2}\nonumber\\ - \frac{2U\left(X'Y''-Y'X''\right)}{{\textit{Rb}}} - \frac{Z''}{{\textit{Fr}}^2}, \end{gather}$$

(3.18)\begin{gather} 0 = \frac{X\left(Y'Z''-Z'Y''\right)+Y\left(Z'X''-X'Z''\right)}{{\textit{Rb}}^2} - \frac{X'Y''-Y'X''}{{\textit{Fr}}^2} + \frac{2UZ''}{{\textit{Rb}}}, \end{gather}

(3.19)\begin{gather} X'^2+Y'^2+Z'^2 = 1. \end{gather}

These equations form a closed system for $X(s)$, $Y(s)$, $Z(s)$, the Cartesian coordinates of the baseline of the jet, and $R(s)$, $U(s)$, which express the local velocity and jet thickness. The factor 3 in front of the second term on the right-hand side of (3.16) is the Trouton number that represents the ratio of extensional to shear viscosity. It appears naturally in the asympotic expansion of the governing equations. In the inviscid limit (${\textit {Re}} \gg 1$), which approximately hold for our experiments, (3.16) can be combined in a single integrand and readily integrated as

(3.20)\begin{equation} U^2+\frac{2U^{1/2}}{{\textit{We}}\sqrt{C_1}}+\frac{2}{{\textit{Fr}}^2}Z-\frac{1}{{\textit{Rb}}^2}\left(X^2+Y^2\right) + C_2 = 0, \end{equation}

where we used $UR^2=C_1$ from (3.15), and $C_1$ and $C_2$ are integration constants. Note that viscous resistances to bending and twisting are also neglected in this case.

The rotation path of the orifice lies on the $XY$ plane, as shown in figure 2(a). We define $X(0)=1$, $Y(0)=0$, $Z(0)=0$, $X'(0)=1$, $Y'(0)=0$ and finally $Z'(0)$ from (3.19), i.e. the jet is issued in the $X$ direction. The nozzle velocity $U_0$ and the radius $R_0$, which have been used for non-dimensionalization, set the base flow boundary conditions as $U(0)=R(0)=1$ and, hence, $C_1=1$. Combining all the prescribed boundary conditions in (3.20) yields $C_2 = 1/{\textit {Rb}}^2 - 2/{\textit {We}} - 1$.

We solve the system (3.15), (3.17)–(3.20), also known as the ‘string’ equations, with the prescribed boundary conditions and the non-dimensional numbers given in table 1 using a fourth-order Runge–Kutta method. The results are summarized in figure 3. The base jet trajectory is in great agreement with the experiments, so it enables the quantification of the projections of the body forces along the trajectory, namely, the centrifugal force, $\,f_{\textit {centrifugal}}(s)$, and gravity force, $\,f_{\textit {gravity}}(s)$, whose expressions are given on the right-hand side of (3.16). Figure 3 depicts the results from the base flow calculations. It is apparent that with the parameters in our experiment, the gravity force is about an order of magnitude smaller than the centrifugal force, which suggests that the torsion along the arc length and distance travelled by the jet in the $Z$ direction are negligibly small, allowing us to study the dynamics from two-dimensional projections onto the $XY$ plane.

Figure 3.

Results from the base flow computation. (a) Comparison of calculated jet trajectory (red line) with a snapshot from an unperturbed jet experiment. (b) The variation of the base state velocity, $U(s)$, and jet radius, $R(s)$, along the jet. (c) The variation of the body forces acting along the jet due to rotation $\,f_{\textit {centrifugal}}=(XX'+YY')/{\textit {Rb}}^2$ and gravity $\,f_{\textit {gravity}}=-Z'/{\textit {Fr}}^2$.

3.3. Nonlinear slender-jet model

To be able to capture the dynamics close to breakup, e.g. the formation of main and satellite droplets, nonlinear simulations are necessary. Once the longitudinal projections of the centrifugal and the gravity forces are accounted for, one can approximate the nonlinear dynamics of the jet along the flow direction with an error of ${O}(\epsilon )$, provided that the full free surface curvature is accounted for (Li et al. Reference Li, Sisoev and Shikhmurzaev2019). In other words, the 3-D spiralling trajectory of the jet can be represented as a quasi-straight jet with locally varying body forces. This also amounts to the implicit assumption that the wavenumber of the perturbations $k$ (scales with $R_0^{-1}$) is much larger than the curvature $\kappa$ and the torsion $\tau$ of the base jet trajectory (both scale with $L^{-1}$) that is valid in the limit $\epsilon \rightarrow 0$. To this end, we implemented a slender-jet approximation, detailed in Eggers & Dupont (Reference Eggers and Dupont1994) and Kamis et al. (Reference Kamis, Eral and Breugem2021), which yields the following unsteady 1-D equations for the jet flow:

(3.21a)$$\begin{gather} \frac{\partial A}{\partial t}+\frac{\partial}{\partial s}\left(u A\right) = 0, \end{gather}$$

(3.21b)$$\begin{gather}\frac{\partial \left(uA\right)}{\partial t} + \frac{\partial}{\partial s}\left(u^2 A\right) = \frac{\partial}{\partial s}\left[A\left(\frac{K}{{\textit{We}}}+\frac{3\epsilon}{{\textit{Re}}}\frac{\partial u}{\partial s}\right)\right] + A\left[\epsilon\left(\,f_{\textit{centrifugal}}+\,f_{\textit{gravity}}\right)\right], \end{gather}$$

(3.21c)$$\begin{gather}K = \frac{2}{\sqrt{4A+A_{\textit{s}}^2}}+\frac{4AA_{\textit{ss}}-2A_s^2} {\left(4A+A_{\textit{s}}^2\right)^{3/2}}. \end{gather}$$

Here $A=R^2$ is the normalized local jet cross-sectional area. Note that the full expression for the curvature is used (i.e. both radii of curvature instead of only the leading order one in the expansion). This is key to the success of the slender-jet framework (Entov & Yarin Reference Entov and Yarin1984; Ambravaneswaran, Wilkes & Basaran Reference Ambravaneswaran, Wilkes and Basaran2002).

Note the inclusion of the $\epsilon$ parameter in front of the centrifugal and gravity force in (3.21b) as they have been taken from (3.16) and are hence normalised with $L$ instead of $R_0$. Equations (3.21) form a closed system for the jet cross-sectional area $A(s,t)$ and jet velocity $u(s,t)$. The initial conditions are implemented as $A(s,0) = R(s)^2$ and $u(s,0) = U(s)$ from the base flow computations detailed in § 3.2. At the nozzle, we input the velocity fluctuations as $u(0,t) = 1+\varepsilon _v \cos {(\omega t)}$ on top of the Gaussian white noise that is mentioned in § 4. A regularized capillary pressure is implemented to carry the simulations beyond the pinch-off point (Driessen & Jeurissen Reference Driessen and Jeurissen2011). A finite difference scheme on a staggered grid is used for solving the system given in (3.21). The time integration is done explicitly by using a three-step Runge–Kutta scheme with an adaptive time stepping. Except for velocities, a central differencing scheme is used for evaluating the variables at their half-step neighbours. For velocities, a higher-order total variation diminishing van Leer scheme is used. More details of the numerical method and validation can be found in Kamis et al. (Reference Kamis, Eral and Breugem2021).

Mechanical vibrations are naturally present in industrial applications of spiralling jets such as in prilling, where the typical size of the rotating perforated bucket is of the order of 50 cm. Even in small-scale laboratory experiments, one can never get fully rid of such vibrations, and they might have a significant effect on the jet breakup. The time series of such noisy perturbations are rather difficult to quantify, however, a good proxy is to represent the vibrations by white noise and to characterize the noise strength in terms of the natural (actuation-free) jet breakup length or breakup time (Pimbley & Lee Reference Pimbley and Lee1977).

To simulate the natural breakup in the slender-jet framework, we add Gaussian white noise to the nozzle velocity as follows:

(3.22)\begin{equation} u(0,t) = 1 + \underbrace{\varepsilon_v \cos{(\omega t)}}_{\text{perturbation}} + \underbrace{W(\Delta t)}_{\text{white noise}}. \end{equation}

The white noise is generated using standard random number generators in Matlab and its standard deviation scales as $S\Delta t$, where $S$ is a strength parameter and $\Delta t$ is the computational time step. This ensures that the forcing due to the white noise is uncorrelated with the time step (Kamis et al. Reference Kamis, Eral and Breugem2021). The strength parameter $S$ is then calibrated using the observed experimental natural (i.e. $\varepsilon _v=0$) breakup length. In the present study we subsequently used $S=300$ in all simulations and the non-dimensional time step $\Delta t$ is adjusted dynamically considering the stability restrictions based on the effects of advection, diffusion and surface tension (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992; Wesseling Reference Wesseling2009).

The comparison of the simulation and experiments for the parameters listed in table 1 and in the presence of only background noise, is summarized in table 2.

Table 2.

Mean and standard deviation of the breakup properties of natural (actuation-free) breakup cases. The breakup length is expressed in terms of the nozzle radius $R_0$ and the drop projected area is represented in terms of $R_0^2$.

4. Comparison of simulations with experiments for jet breakup

Simulations and experiments have been conducted for the parameter settings listed in table 1 and for three different actuation frequencies ($\omega$). The strength of perturbations was varied by varying the velocity perturbation amplitude ($\varepsilon _v$) in the simulations and the oscillation amplitude of the piezo actuator ($\Delta h_{\textit {piezo}}$) in the experiments. The simulation results from the slender-jet model for the jet interface have been projected from the curvilinear to the Cartesian frame using 3.7. Figure 4(a) shows the comparison of a slender-jet simulation and an experiment at $\omega =0.7$ and in figure 4(b) we show the change in the mean intact length with respect to the perturbation amplitude at different frequencies. Note that in figure 4(b) we shifted the axes for $\Delta h_{\textit {piezo}}$ and $\varepsilon _v$ as to give a best match between simulations and experiments. The simulations and experiments agree well in terms of trajectory and the jet shape. The presence of the background noise brings in a threshold for the finite-amplitude perturbations to overcome the noise and control the breakup of the spiralling jet, which can be observed in figure 4(b). The errorbars show the standard deviation around the mean for the time series of the jet intact length, defined as the distance in s from the nozzle exit till the location at which pinch-off takes place. This indicates how much the perturbations are coupled to the system, namely, the less the intact length fluctuates, the more the finite-amplitude perturbations are coupled to the jet, hence better control.

Figure 4.

(a) Comparison of the jet interfaces obtained from simulation (shown in red) and from experiment at a perturbation frequency $\omega =0.7$. The piezoplate displacement is $6 \mathrm {\mu } {\rm m}$ and the velocity amplitude in the simulation is 0.5 %. (b) Mean and standard deviation of the jet intact length as a function of piezoplate displacement in the experiments and velocity perturbation amplitude in the simulations.

To make the comparisons between simulations and experiments conveniently, we show the spiralling jet in a space time plot in figure 5 (without actuation) and figure 6 (with actuation) by stitching consecutive time images together. This also allows us to inspect the motion of the droplets after the pinch-off event, such as capillary oscillations and rotation. The time difference between instants is $\Delta t = 0.8$ ms, during which the orifice rotates ${\rm \pi} /150$ rad. By comparing figures 5 and 6 it is clear that jet actuation helps in achieving a regular jet breakup and reduces the fluctuation in the jet intact length. Also, the comparison of the drop size distributions in figures 5(c) and 6(c) reveal a narrower distribution around the mean projected drop area when the jet is actuated. Overall the agreement between the experiments and simulations is good. Figures 5 and 6 also allow for a qualitative comparison in terms of the late jet dynamics, formation of satellite drops and droplet motion after pinch-off.

Figure 5.

Comparison between experiment and simulation for the natural (actuation-free) breakup case. (a) Experiment. (b) Simulation. (c) Comparison of the drop projected area normalized by the mean.

Figure 6.

Comparison between experiment and simulation in the presence of actuation with $\omega =0.9$. (a) Experiments with a piezoplate displacement of $6\ \mathrm {\mu } {\rm m}$. (b) Simulations with a velocity amplitude 0.5 %. (c) Comparison of the drop projected area normalized by the mean.

For the unperturbed jets, the intact length of the jet fluctuates at a higher standard deviation than their perturbed counterpart. The jets often emit filaments that are a few wavelengths long that fragment into droplets after separating from the main jet. Due to the thinning and stretching of the jet, the high-frequency perturbations that are initially stable and damped, may become unstable downstream. Therefore, a larger window of wavelengths is involved in the instability process as compared with the Rayleigh window for straight jets in the absence of external forces.

In the presence of jet actuation, the jet intact length fluctuations are suppressed and the jet regularly emits a filament that breaks up either at the forward/downstream end or the rear/upstream end. For straight jets, it is known that the location of the breakup within a single filament depends both on perturbation frequency and amplitude (Pimbley & Lee Reference Pimbley and Lee1977). If the breakup occurs at the forward end, the satellite drop merges with the main drop that is following it, within a distance of a few wavelengths. Within a narrow window of amplitudes at a given frequency, the breakup location shifts from the downstream end to the rear end, which is the desired regime for applications such as continuous ink jet printing as the satellite drop merges with the leading main drop (McIlroy & Harlen Reference McIlroy and Harlen2019). The main mechanism of merging is the momentum transfer during the finite time between the two successive pinch-off events, also known as the satellite interaction time (Pimbley & Lee Reference Pimbley and Lee1977). In our simulations and experiments, the typical merging length relative to jet pinch-off in all cases is approximately 2–3 wavelengths. The proportion of the main and satellite droplet volumes are given in figure 7. The volumes are normalized with the jet volume over one perturbation wavelength at the nozzle exit, which is conserved under wavelength stretching downstream due to mass conservation. At a given frequency, the perturbation amplitude determines where the droplets are formed, i.e. lower amplitudes result in droplets formed further downstream. The proportion of the volumes of the main and satellite droplets converge as they are formed further downstream. For a nozzle frequency of $\omega =1.2$, which is initially outside of the Rayleigh window, one can see a larger variation in the drop volume and location. This is attributed to the fact that the unstable frequencies of the white noise gain enough to interfere with the imposed finite-amplitude perturbation. Therefore, the jet volume enclosed by the imposed perturbation wavelength at nozzle exit becomes irrelevant and one can see the droplet volumes are exceeding the volume expected from the imposed perturbation.

Figure 7.

Volumes of the main and satellite droplets obtained from the simulations as a function of their location of first appearance. For each frequency, three points are shown corresponding to different perturbation amplitudes (drops are formed further downstream for smaller amplitudes).

The fact that the droplets stay in the base jet trajectory in the experiments (see figure 6) verifies that it is a good approximation to assume that the jet trajectory is steady in the rotating frame and to separate the calculation of the base flow and the drop formation. After the pinch-off, the droplets travel in a straight line with respect to the stationary frame, preserving their linear momentum.

The satellite interaction time leads to another phenomenon in the context of spiralling jets. With straight jets, the momentum transfer is always in the direction of the flow and the satellite drop is either sped up or slowed down depending on the location of the first pinch-off (forward or rear). When the spiralling jets are captured in the stationary frame (figure 6), one can observe that the capillary oscillation of the droplets is superposed with a rotation around the principal axis. The curved nature of the trajectory results in slight variation of the tangential and the normal vectors from the front end to the rear end of the droplets, which leads to a moment around the principal axis and breaks the symmetry of the detached drop around the tangential direction. The locality of the forces along the flow direction imposes a moment during the satellite interaction time, causing the drop oscillations to be 3-D in nature.

As the perturbations grow and the amplitude becomes finite, the shape deformations of the jet are sufficiently large where the dynamics enter a nonlinear regime. The formation of the main and satellite droplets is a nonlinear phenomenon by nature. The evolution of the velocity perturbations deviates from its linear counterpart (i.e. normal mode), shown in figure 8 for $\omega =0.9$, as the satellite droplets begin to form. It manifests itself as the fluctuations take the form of triangular waves starting from the sinusoidal. On the other hand the nonlinear growth rate is in agreement with the linear estimation until very close to the breakup, which is partly due to the fact that the studied case is close to the inviscid limit, where the time scale for axial momentum diffusion is large compared with the inertio-capillary time scale during the early stage of the growth. In other words, the jet is sufficiently deformed only very close to the breakup for viscosity to play a role.

Figure 8.

Comparison of the linear and nonlinear evolution of the velocity perturbations for $\omega = 0.9$ at different velocity perturbation amplitudes at the nozzle, with decreasing strength from top to the bottom. Here $t_b$ refers to the instant of jet breakup. The dotted lines represent the envelope of the velocity perturbations computed from the linear model detailed in § 5.

5. Analogy with gravity-driven straight jets: approximate base flow solutions, linear stability and self-similarity

The agreement between the experiments and simulations demonstrates the validity of modelling spiralling jets as quasi-straight and slender jets with locally varying pseudo forces provided that epsilon is sufficiently small. These forces, and the variation thereof, can be computed from the base trajectory of the jet as discussed in § 3.2. Once the base flow of the jet and the variation of the body forces have been accounted for, one can use these to study the stability of this flow to perturbations. In this section we will linearize the slender-jet model for a linear spatial stability analysis of the flow to infinitesimally small perturbations imposed at the nozzle exit. We will show that this analysis provides insight in the dynamics underlying jet breakup and in the scaling behaviour for the spatial growth rate and jet breakup.

To study how perturbations evolve along the jet, one needs to pay attention to the disparate length scales present in this problem. In § 3.2 where the base flow and the trajectory was calculated in the limit $\epsilon \rightarrow 0$, all lengths are expressed in terms of the rotating arm $L$ and the base flow equations represent the jet as a ‘string’ whose thickness is $\textit {O}(\epsilon L)$. When we consider capillary waves along the jet, their wavelengths are of the same order as the jet radius so their growth should be expressed along $s/(\epsilon L)$, i.e. in units of jet radius instead of the rotating arm. Similarly, while the advective scale of the base flow is $t_{\textit {adv}}=L/U_0$, the relevant advective time for the perturbations is $R_0/U_0 = t_{\textit {adv}}/\epsilon$. In what follows, we will express all the lengths in terms of the nozzle radius and the length of the rotating arm $L$ will only appear in the expression of the body forces. Furthermore, we introduce an effective gravity $g_e\sim \textit {O}(\varOmega ^2L+g)$ that measures the combined effect of rotational and gravitational forces. Using this, we define the effective capillary length as $l_c = \sqrt {\gamma /\rho g_e}$ and an effective Bond number ${\textit {Bo}}_e=(R_0/l_c)^2=\rho g_e R_0^2/\gamma$. In terms of (We, Rb, Fr) given in table 1, the latter can be expressed as

(5.1)\begin{equation} {\textit{Bo}}_e({\textit{We}},{\textit{Fr}},{\textit{Rb}}) = {\textit{We}}\left(\frac{\epsilon}{{\textit{Rb}}^2}+\frac{\epsilon}{{\textit{Fr}}^2}\right). \end{equation}

Using (5.1), one could follow the same line of gain analysis developed for straight jets in Le Dizès & Villermaux (Reference Le Dizès and Villermaux2017). For the given experimental case with the non-dimensional numbers presented in table 1, the variation of the body forces along the jet is minor; see figure 3(c). The centrifugal forces are very close to their nozzle value throughout the entire jet and gravitational forces, which would make the jet travel in the $Z$ direction and cause torsion, are an order of magnitude smaller than the centrifugal forces. In this case, the spiralling jet thus behaves like a gravity-driven straight jet with the gravitational acceleration replaced by the initial centrifugal acceleration at the nozzle exit, namely, $U(dU/ds) = {\textit {Bo}}_e/{\textit {We}}$, as the viscous and curvature effects decay rapidly after the nozzle exit. One thus has

(5.2a)$$\begin{gather} U(s) = \sqrt{\frac{2{\textit{Bo}}_e}{{\textit{We}}}s+1}, \end{gather}$$

(5.2b)$$\begin{gather}R(s) = \left(\frac{2{\textit{Bo}}_e}{{\textit{We}}}s+1\right)^{{-}1/4}. \end{gather}$$

A comparison between the ‘free-fall’ solution (5.2) and the solution of the base flow equations (3.20) is given in figure 9.

Figure 9.

The comparison between the base flow solutions from (3.15), (3.17)–(3.20) and the approximate ‘free-fall’ solutions given in (5.2).

This assumption can be made a priori if ${\textit {We}}\gg {\textit {Bo}}_e$ for ${\textit {We}}$ large enough to stay in the jetting regime (Clanet & Lasheras Reference Clanet and Lasheras1999). The arc length $s$ can be recast as a slowly varying dimension

(5.3a,b)\begin{equation} \xi = \frac{s}{s_0}+1, \quad s_0\,{\rm d}\xi = {\rm d}s, \end{equation}

for large $s_0$, where $s_0 = {\textit {We}}/2{\textit {Bo}}_e$ is the variation scale (Le Dizès & Villermaux Reference Le Dizès and Villermaux2017; Eggers & Villermaux Reference Eggers and Villermaux2008), such that (5.2) becomes

(5.4a)$$\begin{gather} U(\xi) = \sqrt{\xi}, \end{gather}$$

(5.4b)$$\begin{gather}R(\xi) = \xi^{{-}1/4}, \end{gather}$$

respectively.

We now consider perturbations to the base flow of the following form:

(5.5)\begin{equation} (\tilde{u},\tilde{R}) = (u(\xi),r(\xi)) \exp\left({{\rm i}\left(s_0\int_1^\xi k(\xi')\,{\rm d}\xi'- \omega t\right)}\right). \end{equation}

Here $k = k_r + \text {i}k_i$ is the complex wavenumber, $k_r=2{\rm \pi} /\lambda$ is the real wavenumber with $\lambda$ the perturbation wavelength (normalised with $R_0$), $-k_i$ is the spatial growth rate and $\omega$ is the real perturbation frequency imposed at the nozzle. Plugging $U+\tilde {u}$ and $R+\tilde {R}$ into (3.21a)–(3.21c) and considering the terms linear in $(\tilde {u},\tilde {R})$ up to $\textit {O}(1/s_0)$ yields

(5.6a)$$\begin{gather} \left(kU-\omega\right)r(\xi) + \frac{kR}{2}u(\xi) = \textit{O}\left(\frac{1}{s_0}\right), \end{gather}$$

(5.6b)$$\begin{gather}\left(kU-\omega-\frac{3\text{i}k^2}{{\textit{Re}}}\right)u(\xi) + \frac{k}{R^2{\textit{We}}}\left(k^2 R^2 -1\right)r(\xi) = \textit{O}\left(\frac{1}{s_0}\right), \end{gather}$$

which is analogous to equation (2.16) in Le Dizès & Villermaux (Reference Le Dizès and Villermaux2017). This approach is also known as the WKBJ (Wentzel–Kramers–Brillouin–Jeffreys) method (Bender & Orszag Reference Bender and Orszag1978), a perturbation method using expansions over the small parameter $1/s_0$ in this case. This approach is well utilized and tested in Javadi et al. (Reference Javadi, Eggers, Bonn, Habibi and Ribe2013) and Le Dizès & Villermaux (Reference Le Dizès and Villermaux2017). The variations of the amplitudes $(u(\xi ),r(\xi ))$ and the base states only show up at higher orders, which means for large $s_0$ (i.e. ${\textit {We}}\gg {\textit {Bo}}_e$), the length scale of base flow variations are small compared with the perturbation wavelength. The dispersion relation is found from the non-trivial solution of (5.6),

(5.7)\begin{equation} \frac{k^2}{2R {\textit{We}}}\left(k^2 R^2 -1\right) - \left(kU-\omega\right)^2 + \frac{3\text{i}k^2}{{\textit{Re}}}\left(kU-\omega\right) = 0. \end{equation}

This equation governs the convective instability of a jet to infinitesimally small perturbations at the nozzle exit at a given frequency $\omega$. Note that the dispersion relation (5.7) is identical to that for a straight jet in the absence of body forces, except that here $U$ and $R$ are the local values instead of the nozzle values.

The solution to the dispersion relation given in (5.7) is shown in figure 10. The local wavenumber, $k_r$, is the real part of the solution to the complex equation and it is decreasing from its nozzle value due to the thinning of the jet, also known as wavelength stretching (Eggers & Villermaux Reference Eggers and Villermaux2008). Based on the nozzle values, the frequency that has the maximum growth rate at the nozzle is $\omega _{\textit {ini}}\approx 0.718$. But as the body forces stretch the perturbed filaments, their growth rate also changes locally which leads to different gains along the jet, shown in figure 11(a). So a net gain approach is preferable to study the most unstable frequencies (Javadi et al. Reference Javadi, Eggers, Bonn, Habibi and Ribe2013; Le Dizès & Villermaux Reference Le Dizès and Villermaux2017). The local gain along the jet is then expressed as

(5.8a,b)\begin{equation} G(\xi) = e^{S(\xi)}, \quad S(\xi) = s_0\int_1^\xi -k_i(\xi')\,{\rm d}\xi', \end{equation}

where $s_0 = {\textit {We}}/2{\textit {Bo}}_e \sim 0.5\lvert \,{\rm d}U/{\rm d}s\rvert _0^{-1}/(R_0/U_0)$ can be interpreted as the ratio of the initial stretching time to the initial advective time.

Figure 10.

The local wavenumber ($k_r R_0$) and spatial growth rate ($-k_i R_0$) along the jet. The dotted line indicates the upper boundary of the local Rayleigh window of unstable perturbations.

Figure 11.

(a) Integrated gain for different frequencies of perturbation at the nozzle. The dashed vertical line indicates the location $\xi =2.35$ that corresponds to the location of the mean jet intact length obtained from the experiments in the absence of imposed perturbations, also given in table 2. (b) The location to reach the gain necessary for breakup versus the frequency of perturbation at the nozzle.

Figure 11(a) shows the gain along the jet at different frequencies and for infinitesimal amplitude, based on (5.8a,b). The disturbances start to gain only when their local wavelength enters the Rayleigh window, which is indicated by the dotted line in figure 10. Quantification of the infinitesimal disturbances is a difficult practice both experimentally and numerically. When the gain reaches a sufficiently high value, which is assumed to be $G\approx e^7$ in Le Dizès & Villermaux (Reference Le Dizès and Villermaux2017), the transition to breakup occurs. The location of the transition $\xi _t$ for different values of threshold gain is given in figure 11(b). For our experimental case with the parameters given in table 1, the assumed transition gain of $S_t\approx 7$ agrees well with the average natural breakup length for the unperturbed case of our experiments that is around 95 times the radius of the jet; see table 2. For the choice of $S_t\approx 7$, the frequency that gains the fastest is around $\omega \approx 0.9$.

Furthermore, self-similarity can be deduced from the dispersion relation. We will first decompose the solution of the complex equation (5.7) as

(5.9a)$$\begin{gather} k_rR(\xi) = \frac{\omega R}{U}= \omega\xi^{{-}3/4}, \end{gather}$$

(5.9b)$$\begin{gather}-k_iR(\xi) = q\left(k_rR,{\textit{Oh}}\right){\textit{We}}_l^{{-}1/2} = q\left(\omega \xi^{{-}3/4},{\textit{Oh}}\right){\textit{We}}^{{-}1/2}\xi^{{-}3/8}, \end{gather}$$

where $q$ is a non-dimensional function of the local wavenumber, the perturbation frequency and the Ohnesorge number, and ${\textit {We}}_l = We\xi ^{3/4}$ is the local Weber number obtained by using the base flow solutions (5.4). Looking into the convective behaviour of the perturbations of $\omega >1$, one can see that the initial wavelength of these perturbations are shorter than the lower boundary of the Rayleigh window (i.e. $k_r(\xi )R\leq 1$ as shown with a dashed curve in figure 10), so these perturbations do not gain (or even damped with higher viscosity) until the wavelength is stretched enough to enter the Rayleigh window. This motivates us to do another rescaling on the flow direction, $\xi$, as

(5.10)\begin{equation} \theta = \xi\omega^{{-}4/3}, \end{equation}

such that $k_rR(\theta )=\theta$ and instabilities grow for $\theta >1$. From here, we can express (5.9) explicitly in terms of $\theta$ using the base flow solutions (5.4),

(5.11a)$$\begin{gather} k_r\omega^{{-}1/3} = \theta^{{-}1/2}, \end{gather}$$

(5.11b)$$\begin{gather}-k_i\omega^{1/6}{\textit{We}}^{1/2} = q\left(\theta^{{-}3/4},{\textit{Oh}}\right)\theta^{{-}1/8}=p(\theta,{\textit{Oh}}), \end{gather}$$

where $k_r$ and $k_i$ are now scaled with parameters defined at the nozzle exit. Using these rescaled solutions, we can collapse the curves given in figure 10 onto a single curve, which is shown in figure 12(a). We can also write the integral for the gain $S(\xi )$ given in (5.8a,b) in terms of an integral over the similarity variable $\theta$ using (5.11b),

(5.12)\begin{equation} S(\theta) = \frac{s_0\omega^{7/6}}{\sqrt{{\textit{We}}}}\int_1^\theta p(\theta',{\textit{Oh}})\,{\rm d}\theta'. \end{equation}

Figure 12.

(a) Rescaled curves of the solution to the dispersion relation given in (5.7). (b) The gain associated with the rescaled curve.

This is plotted in figure 12(b) for different values of $\omega$. For a threshold gain of $S_t\approx 7$, we can then deduce the rescaled threshold distance for jet breakup, $\theta _{t}$, numerically from (5.12). This is shown in figure 13(a). In figure 13(b) the threshold distance, $\xi _t = \theta _t\omega ^{4/3}$, and the wavelength at breakup calculated as $\lambda (\xi _t) = ({2{\rm \pi} }/{\omega })\xi _t^{1/2}$ are shown. The prefactor $s_0/\sqrt {{\textit {We}}}\sim 0.5\lvert \, {\rm d}U/{\rm d}s\rvert _0^{-1}/\sqrt {\rho R_0^3/\gamma }$ in (5.12) can be interpreted as the ratio of the initial stretching time to the initial inertio-capillary time.

Figure 13.

(a) The rescaled distance for transition to breakup, $\theta _t$, based on (5.12) as a function of frequency. (b) Transition distance, $\xi _t$, based on (5.8a,b) and the wavelength at the transition distance as a function of frequency calculated as $\lambda = 2{\rm \pi} /k_r(\xi _t)$.

As nonlinear effects present in the experiments are a priori excluded from the linear stability analysis, a direct quantitative comparison between the linear stability predictions and the experiments or the nonlinear slender-jet model would be flawed. Nevertheless, we made a quantitative comparison of the linear and nonlinear predictions for the evolution of the perturbations for the non-dimensional numbers given in table 1. To that end, we performed a small study where we compare the results from linear stability analysis with the nonlinear slender-jet model for two different perturbation amplitudes and a range of frequencies.

Results for the jet breakup length are given in figure 14 where the breakup length $L_b$ is normalized by the minimum breakup length, $L_{b,min}$. One can see that although the value of the most unstable wavelength does agree between the linear and nonlinear predictions, the behaviour at other wavelengths differs significantly. At high frequencies (i.e. short wavelengths) and $\varepsilon _v = 1\,\%$, the nonlinear prediction for the jet breakup length does follow the $\omega ^{4/3}$ scaling, though slightly underestimating the breakup length. However, at low frequencies (i.e. long wavelengths) there is a significant overestimation of the breakup length by the linear analysis, as for $\omega \rightarrow 0$, the perturbation wavelengths become eventually large enough to be comparable to the streamwise variation scale of the base flow. To correct for this, it would be necessary to include higher-order terms in (5.6), starting with terms up to $O(1/s_0^2)$ and increasing the order for smaller $\omega$. In our study the focus is on the frequencies that are in the vicinity of the most unstable frequencies, corresponding to the larger frequencies and smaller perturbation wavelengths for which the local plane-wave approximation underlying (5.6) is valid.

Figure 14.

Comparison of the breakup length from linear stability theory and from nonlinear simulations for two different velocity fluctuation amplitudes, $\varepsilon _v$, as a function of the imposed perturbation frequency. To ease the comparison, the breakup length has been normalised with the minimum breakup length at the most unstable frequency. Values for $L_{b,min}$ are 91$R_0$ from the linear stability, and 90.6$R_0$ and 52.7$R_0$ from the nonlinear simulations for $\varepsilon _v = 0.1\,\%$ and $\varepsilon _v = 1\,\%$, respectively.