1. Introduction
The breakup of liquid jets into droplets widely occurs in natural phenomena and industrial applications, such as water faucets and fountains (Eggers & Villermaux Reference Eggers and Villermaux2008), archerfish predation (Gerullis & Schuster Reference Gerullis and Schuster2014), inkjet printing (Lohse Reference Lohse2022), and fuel atomisation in machines (Ren et al. Reference Ren, Wang, Xiang, Zhao and Zheng2019). In these scenarios, the liquid phase is injected from a round nozzle into a static air environment, forming a cylindrical jet and eventually disintegrating into a chain of droplets due to interfacial instability. It should be emphasised that the morphology of droplets is of vital importance to the efficiency of applications. For example, in the case of tin droplets generated in extreme ultraviolet lithography, the size, homogeneity and separation distance of droplets must be precisely controlled to ensure the highly repeatable stimulation of tin droplets under strong laser pulses (Sun et al. Reference Sun, Qi, Wang and Zuo2024). For the fabrication of microcapsules in pharmacy, the droplet productivity and size homogeneity are important for process optimisation and functional consistency (Gañán-Calvo et al. Reference Gañán-Calvo, Montanero, Martin-Banderas and Flores-Mosquera2013). Therefore, controlling the interface breakup of liquid jets is of great significance, as it facilitates the active generation of droplets.
The investigation of liquid jet breakup dates back to Savart (Reference Savart1833), who conducted an experimental study and found that the development of interface perturbations led to the eventual disintegration of the liquid jet. The growth rate of perturbations could be measured based on the temporal development of the perturbation wave crest through stroboscopic photography (Donnelly & Glaberson Reference Donnelly and Glaberson1966). Satellite droplets have also been observed to accompany the main droplets, which were caused by the strong nonlinear characteristics of the jet interface during breakup (Hoeve et al. Reference Hoeve, Gekle, Snoeijer, Versluis, Brenner and Lohse2010). The occurrence of satellite droplets was found to be closely related to the liquid jet viscosity. As the jet velocity increases, shear stress could play a more significant role, leading to a transition from Rayleigh breakup (dominated by surface tension) to the first wind-induced regime, the second wind-induced regime and finally to the atomisation regime (Lin & Reitz Reference Lin and Reitz1998; Zhan et al. Reference Zhan, Kuwata, Maruyama, Okawa, Enokia, Aoyagib and Takata2020). At relatively low liquid velocity, the interface was observed to break up right at the nozzle exit, corresponding to the dripping mode (Vihinen, Honohan & Lin Reference Vihinen, Honohan and Lin1997; Clanet & Lasheras Reference Clanet and Lasheras1999). Numerical simulations have also been conducted to reveal the dynamics of jet breakup, providing more quantitative details that are difficult to measure in experiments (Subramani et al. Reference Subramani, Yeoh, Suryo, Xu, Ambravaneswaran and Basarana2006; Zhou, Yue & Fenga Reference Zhou, Yue and Fenga2006; Huang, Bao & Qian Reference Huang, Bao and Qian2023). Linear instability analysis is a powerful tool for studying the perturbation development of liquid jets from the theoretical perspective (Lin Reference Lin2003). The pioneering work of Plateau (Reference Plateau1873) and Rayleigh (Reference Rayleigh1878) indicated that the surface tension force could lead to the axisymmetric instability of an inviscid jet. The liquid viscosity was found to weaken jet instability and modulate the wavelength of the most unstable perturbation (Tomotika Reference Tomotika1935; Chandrasekhar Reference Chandrasekhar1961). The non-axisymmetric helical mode of jet instability tends to occur with the addition of external fields, such as the external driving gas stream (Si et al. Reference Si, Li, Yin and Yin2009, Reference Si, Li, Yin and Yin2010), the electric field (Loscertales et al. Reference Loscertales, Barrero, Guerrero, Cortijo, Marquez and Gañán-Calvo2002; Li, Yin & Yin Reference Li, Yin and Yin2009) and the circumferential angular velocity (Gallaire & Chomaz Reference Gallaire and Chomaz2003; Kubitschek & Weidman Reference Kubitschek and Weidman2007; Xu et al. Reference Xu, Mu, Xin, Qiao, Zhao and Si2024). To date, several review articles have been published summarising the instability characteristics of liquid jets under various conditions, both from physical and theoretical perspectives (Lin & Reitz Reference Lin and Reitz1998; Eggers & Villermaux Reference Eggers and Villermaux2008; Montanero & Gañán-Calvo Reference Montanero and Gañán-Calvo2020).
Due to random perturbations on the jet interface, the breakup of liquid jets into droplets always presents relatively poor repeatability. For example, the droplet size usually falls within some ranges and the satellite droplets always form accompanying the main droplet. Moreover, the generating frequency and average size of droplets are hard to be adjusted actively. These shortcomings limit the efficiency of practical applications, highlighting the need for on-demand control of jet breakup and droplet production. Currently, multiple methods have been introduced to modulate the breakup of liquid jets, such as the addition of laser illumination (either continuous laser (Liu et al. Reference Liu, Wang, Gao, Huang, Tang, Zhao and Deng2021) or periodic pulsed laser (Zhao et al. Reference Zhao, Wan, Chen, Chao and Xu2021)), ultrasonic field (Kudryashova et al. Reference Kudryashova, Shalunov, Titov, Dorovskikh and Nesterov2023), thermal field (Basaran, Gao & Bhat Reference Basaran, Gao and Bhat2013; Kamis, Eral & Breugem Reference Kamis, Eral and Breugem2021) and mechanical vibration (She et al. Reference She, Fang, Hu, Su and Fu2022; Luo et al. Reference Luo, Lyu, Qi and Li2023). Among these methods, the mechanical vibration acts as the simplest one and has been most widely used. The mechanical vibration is able to provide an additional perturbation to the liquid jet, which competes with the jet intrinsic perturbation and is able to modulate the jet breakup behaviour under certain conditions. Moallemi, Li & Mehravaran (Reference Moallemi, Li and Mehravaran2016) have found in experiments that a train of droplets with uniform size can be generated by applying mechanical vibration with a sinusoidal waveform, and the droplet size and jet breakup length could be modulated within a certain range by changing the vibrating frequency and amplitude, respectively. These findings have also been comfirmed by their numerical simulations with adaptive mesh refinement. The frequency range for regular jet breakup with uniform droplets was experimentally determined by Rohani, Jabbari & Dunn-Rankin (Reference Rohani, Jabbari and Dunn-Rankin2010), considering the variation of jet velocity. It should be noted that, apart from free liquid jets in static ambient gas, the mechanical vibration can also modulate the breakup of liquid jets in complex flow situations, such as liquid jets driven by electric fields (Duan et al. Reference Duan, Yang, Li, Lojewski and Deng2013), coflowing gas flows (Xu et al. Reference Xu, Zhu, Mu, Huang and Si2022) and immiscible liquid streams (Zhu, Tang & Wang Reference Zhu, Tang and Wang2016; Yang et al. Reference Yang, Qiao, Mu, Zhu, Xu and Si2019). Liquid jets were also found to generate periodic wrinkled structures instead of droplets under very low interfacial tension between the jet and the surrounding liquid (Sauret, Spandagos & Shum Reference Sauret, Spandagos and Shum2012; Sauret & Shum Reference Sauret and Shum2012).
For the liquid jet breakup in a static gas environment, it has been found that the combination of multiple sinusoidal perturbations with different wavelengths can lead to more complex behaviours of droplet generation and subsequent evolution. It was observed by Orme & Muntz (Reference Orme and Muntz1990) that the superposition of two sinusoidal perturbations could lead to the initial jet breakup into droplets at distances commensurate with the faster frequency of the disturbance, and droplets from two adjacent half-periods of the modulation cycle would eventually merge to form one large drop as they travel downstream. The merging process of multiple droplets was closely related to the relative amplitude between the two disturbances (Rohani et al. Reference Rohani, Jabbari and Dunn-Rankin2010). Driessen et al. (Reference Driessen, Sleutel, Dijksman, Jeurissen and Lohse2014) applied the superposition of two perturbations with similar growth rates, where the wavelengths of the perturbations were also close to the wavelength of the most unstable mode. It was found that the liquid jet first degenerated into multiple droplets, which then merged to form final uniform droplets with large separation distances. The period of the final uniform droplets corresponds to the shortest common period of the two perturbations.
Despite the existing studies on the forced breakup of liquid jets, a systematic study of the jet response characteristics over a wide range of perturbation frequencies, amplitudes and jet velocities is still desired, and the transition criteria among different response modes remain to be explored. Moreover, beyond the regular breakup region, the merging dynamics of inhomogeneous droplets with different sizes and the methods for regulating droplet merging also motivate our interest. In this work, we focus on the forced breakup of liquid jets in ambient gas surroundings through direct numerical simulations, in which a velocity disturbance is applied at the jet inlet to modulate the breakup. Theoretical analysis based on the linear instability model is also employed to quantitatively analyse the numerical results. The paper is organised as follows. In § 2, we present our numerical methods based on the diffuse interface method. Section 3 demonstrates the response dynamics of the liquid jet with varying actuation frequency, amplitude and jet velocity. The development of interface perturbations under different response modes is also analysed in this section. The theoretical analysis is provided in § 4. Section 5 investigates the droplet merging dynamics under the mode where multiple droplets are generated in one period, and the manipulation strategy for this mode by adding a pulse to the sinusoidal perturbation wave is studied in § 6. The main conclusions of this paper are presented in § 7.
2. Numerical methodology
The evolution of a liquid jet with a periodically modulated flow rate at the inlet is considered in our study, as sketched in figure 1. The boundary conditions for the calculation are also given in the figure. In the numerical simulations, a liquid jet (with density
$\rho _1$
, dynamic viscosity
$\mu _1$
, initial radius
$R$
, flow rate
$\bar Q_1(\bar t)$
and inlet velocity
$\bar U_1(\bar t)$
, where
$\bar t$
denotes time) flows from the left-side inlet into the static gas surroundings (with density
$\rho _2$
and dynamic viscosity
$\mu _2$
). As the liquid jet evolves downstream, it eventually breaks up into droplets due to the interfacial tension force. Our study mainly focuses on the situation where a sinusoidal waveform is imposed on the inlet velocity of the jet, as illustrated by the inset of figure 1. The detailed description of the inlet velocity is
\begin{eqnarray} \bar U_1(\bar t)=\bar U_1 \big(1+A \sin (2{\pi }\bar f \bar t)\big), \end{eqnarray}
where
$\bar U_1$
is the average velocity of the jet, and
$\bar f$
and
$A$
represent the frequency and amplitude of the perturbations, respectively.
Figure 1. Sketch of the axisymmetric computational domain for numerical simulations, where
$\bar U_1(\bar t)$
represents the oscillating inlet velocity of the liquid jet. The inset sketches the temporal evolution of inlet velocity, where
$\bar U_1$
is the jet average velocity, and
$A$
and
$\bar f$
represent the amplitude and frequency of perturbations, respectively.
In numerical simulations, the diffuse interface method is used to distinguish between the liquid and gas phases. In this method, an interface with finite thickness is assumed. The interface between the liquid and gas is represented by a volume fraction
$C$
, which varies continuously from 1 to 0 across the interface. The finite thickness characteristic of the liquid interface provides a smooth variation of stress, avoiding the stress singularity and thus providing a more accurate calculation of the interfacial tension force compared with sharp interface methods, such as level-set or volume-of-fluid (Ding, Spelt & Shu Reference Ding, Spelt and Shu2007). The interface evolution is governed by the Cahn–Hilliard equation (Jacqmin Reference Jacqmin1999):
\begin{equation} \frac {\partial C} {\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol {u} C) = \frac {1}{\textit{Pe}}{\nabla} ^2 \psi , \end{equation}
where
$\textit{Pe}$
is the Péclet number and the chemical potential
$\psi$
is defined as
\begin{equation} \psi = C^{3} - 1.5 C^{2}+ 0.5 C - Cn ^{2} {\nabla} ^2 C, \end{equation}
where the Cahn number
$Cn$
is a direct measurement of the thickness of the diffuse interface. Generally,
$Cn \sim \Delta x$
to better resolve the interface and ensure accurate surface tension calculation (Ding et al. Reference Ding, Spelt and Shu2007), where
$\Delta x$
is the mesh size. A larger value of
$Cn$
results in a thicker diffuse interface but a more accurate calculation of surface tension force, as more meshes are used to resolve the interface. Based on our previous numerical experience (Mu et al. Reference Mu, Qiao, Si, Chen and Ding2021, Reference Mu, Zhang, Si and Ding2022, Reference Mu, Qiao, Ding and Si2023), we set
$Cn=0.5\Delta x$
. The Péclet number
$\textit{Pe}$
represents the relative significance of convective fluxes compared with diffusive fluxes. It has been suggested by Magaletti et al. (Reference Magaletti, Francesco, Chinappi, Marino and Casciola2013) that the diffuse interface can approach the sharp interface limit as
$Cn$
vanishes, with
$\textit{Pe} \sim Cn^{-1}$
. Therefore, we choose
$\textit{Pe} = 1/Cn$
in our simulations.
The motion of the liquid jet and gas is governed by the dimensionless Navier–Stokes equations:
\begin{align} \rho \left(\frac {\partial \boldsymbol {u}}{\partial t}+\boldsymbol {u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol {u} \right) = -\boldsymbol{\nabla }p +\frac {1}{\textit{Re}}\boldsymbol{\nabla }\boldsymbol{\cdot }\big[\mu \big(\boldsymbol{\nabla }\boldsymbol {u}+\boldsymbol{\nabla }\boldsymbol {u}^T \big) \big]+\frac {\boldsymbol{f}_{\kern-1.5pt s}}{\textit{We}}, \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {u} = 0, \\[9pt] \nonumber \end{align}
where
$\boldsymbol{f}_{\kern-1.5pt s}=6\sqrt {2}\psi \boldsymbol{\nabla }C / Cn$
denotes the surface tension force, and
$\rho =C + (1- C)r_{\rho }$
and
$\mu = C + (1- C)r_{\mu }$
are the dimensionless averaged density and viscosity, in which
$r_{\rho }=\rho _2/\rho _1$
and
$r_{\mu }=\mu _2/\mu _1$
are the density and viscosity ratios between the gas phase and the liquid phase, respectively. In the liquid phase (i.e.
$C=1$
),
$\rho =1$
and
$\mu =1$
, while in the gas phase (
$C=0$
),
$\rho =r_{\rho }$
and
$\mu =r_{\mu }$
. At the interface, the values of
$\rho$
and
$\mu$
change gradually with the variation of
$C$
. Taking the liquid jet as the characteristic phase,
$R$
as the characteristic length and
$\bar U_1$
as the characteristic velocity, the Reynolds and Weber numbers in (2.4) can be defined as
$\textit{Re}=\rho _1 \bar U_1 R/\mu _1$
and
$\textit{We}=\rho _1 \bar U_1^2 R/\sigma$
. As the variation of
$\bar U_1$
will change the values of
$\textit{Re}$
and
$\textit{We}$
simultaneously, we also define the Ohnesorge number as
$\textit{Oh}=\mu _1/\sqrt {\rho _1 \sigma R}$
for the convenience of analysis. The
$\textit{Oh}$
can also be written as
$\textit{Oh}=\sqrt {\textit{We}}/\textit{Re}$
. The characteristic time and frequency correspond to
$R/\bar U_1$
and
$\bar U_1/R$
, respectively. Therefore, the dimensionless form of the inlet velocity
$U_1(t)$
can be written as
\begin{eqnarray} U_1(t)=1+A \sin (2{\pi }ft), \end{eqnarray}
where
$t$
denotes the dimensionless time, and
$A$
and
$f$
represent the dimensionless amplitude and frequency of perturbations, respectively.
The numerical simulations are conducted in axisymmetric cylindrical coordinates (denoted by
$r{-}z{-}\theta$
) with a uniform Cartesian mesh, where
$r$
,
$z$
and
$\theta$
denote the radial, axial and azimuthal directions, respectively. Due to the axisymmetric geometry of a liquid jet, the velocity component at
$\theta$
-direction equals to zero invariably, and only the jet evolutions in the
$r$
- and
$z$
-directions need to be calculated. The computational length in the
$r$
-direction is set to
$8R$
and in the
$z$
-direction, it varies from
$100R$
to
$300R$
, depending on the perturbation amplitude and frequency, which significantly affect the jet breakup length. The interface is represented by the contour of
$C=0.5$
. The mesh size used in the calculation is
$\Delta x=0.05R$
, unless otherwise stated. Based on our previous study (Mu et al. Reference Mu, Si, Li, Xu and Ding2018), this mesh size is fine enough to resolve the liquid jet and droplets. The boundary conditions for the flow velocities are as follows:
$v=0, \partial {u}/\partial {r}=0$
at the axial-symmetric axis of
$r=0$
, where
$u$
and
$v$
are the flow velocities in the
$z$
- and
$r$
-directions, respectively;
$\partial {u}/\partial {z}=0$
and
$\partial {v}/\partial {t}+v \boldsymbol{\cdot }(\partial {v}/\partial {z})=0$
at the right-side outlet;
$\partial {v}/\partial {r}=0$
and
$\partial {u}/\partial {t}+u \boldsymbol{\cdot }(\partial {u}/\partial {r})=0$
at the upper-side outlet;
$u=U_1(t)$
and
$v=0$
at the left-side inlet, where the specific form of
$U_1(t)$
is given by (2.6).
Figure 2. Validation of numerical simulations with experiments of Moallemi et al. (Reference Moallemi, Li and Mehravaran2016). (a) Experimental and numerical interface profiles of an actuated water jet at fixed values of
$A=1 \times 10^{-4}, f=0.111, \textit{Re}=367, \textit{Oh}=0.0105, r_{\rho }=0.001$
and
$r_{\mu }=0.025$
, respectively, where different mesh sizes are considered. (b) Comparison of jet breakup length
$L_j$
(under fixed mesh size of
$\Delta x=0.05R$
) with varying perturbation amplitude of
$A$
.
The comparison of liquid jet profiles between numerical simulations and experimental results is given in figure 2(a). The experiments were conducted by Moallemi et al. (Reference Moallemi, Li and Mehravaran2016), where a free water jet with a radius of
$R$
= 0.125 mm and a mean velocity of
$\bar U_1-2.398\,\text{m s}^{-1}$
was injected into a quiescent air environment. The imposed perturbation frequency was 2609 Hz and the perturbation amplitude was
$10^{-4} \bar U_1$
. The liquid jet breakup presented an axisymmetric manner in experiments, indicating the rationality of axisymmetric numerical simulation. The dimensionless parameters in numerical simulations correspond to
$\textit{Re}=367$
,
$\textit{Oh}=0.0105$
,
$r_{\rho }=0.001$
,
$r_{\mu }=0.025$
,
$f=0.111$
and
$A=10^{-4}$
. To present the mesh convergence of the calculation, different mesh sizes (
$\Delta x=0.05R$
,
$0.025R$
and
$0.02R$
) are employed. Clearly, a reasonable agreement can be obtained with respect to the jet length, and the position and shape of the generated droplets. Moreover, as the mesh size is finer than
$0.05R$
, there is little influence of mesh size on jet breakup. In this work, to provide more quantitative validation, we further vary the actuation amplitude
$A$
over a wide range and compare the dimensionless jet breakup length (denoted as
$L_j$
) from our simulations with the experimental and numerical results of Moallemi et al. (Reference Moallemi, Li and Mehravaran2016). The results are presented in figure 2(b). In our simulations, the jet breakup length is obtained by tracking the temporal evolution of the jet tip position, and the details will be provided in figure 4. It is notable that as the liquid jet breaks up into droplets periodically, the jet breakup length remains almost constant under fixed values of
$A$
and
$f$
, resulting in almostly negligible statistical error bars, as shown by Moallemi et al. (Reference Moallemi, Li and Mehravaran2016). For the convenience of data comparison, the statistical error bars are omitted in figure 2. Clearly, a good agreement can be reached between our simulations and the results of Moallemi et al. (Reference Moallemi, Li and Mehravaran2016) as
$A$
varies. Based on the current values of jet radius and velocity (non-dimensionalised by fixed
$\textit{Re}=367$
,
$\textit{Oh}=0.0105$
,
$r_{\rho }=0.001$
and
$r_{\mu }=0.025$
, regarded as the reference state), this study would systematically evaluate the influence of perturbation frequency
$f$
and amplitude
$A$
on the dynamics of liquid jet breakup and droplet generation, as shown by §§ 3.1 and 3.2. The effect of velocity on jet breakup would also be analysed since the variation of jet velocity could change the natural breakup characteristics of the unactuated jet, thus modulating the competition mechanism between the external actuation and the intrinsic perturbation. As the jet velocity changes comparing with the reference state, the numerical calculation would modulate the dimensionless jet velocity at fixed
$\textit{Re}$
instead of varying
$\textit{Re}$
under constant dimensionless jet velocity of unity. The convenience of this approach will be elucidated in detail in § 3.3.
3. Dynamics of forced jet breakup
3.1. Response modes and phase diagram
As the external perturbation is applied to the liquid jet, the jet breakup behaviour is closely related to the perturbation amplitude
$A$
and frequency
$f$
. Figure 3(a) shows four typical response modes of jet breakup as
$f$
and
$A$
vary, and the dynamic evolutions of these modes are also given as supplementary movie 1 is available at https://doi.org/10.1017/jfm.2025.10833. The phase diagram of different modes is further depicted in figure 3(b). For the first response mode, the jet breaks up randomly and the size of the droplets falls within a certain range. This mode is defined as the irregular breakup mode (referred to as the I mode hereafter). It should be noted that for the free liquid jet without external actuation, the jet breakup also exhibits the I mode due to the random perturbations on the jet interface. The second type of response mode leads to droplets generation with multiple sizes in a periodic manner, with the period for the droplet series equal to
$1/f$
. This is referred to as the multiple breakup mode (M mode). In the third response mode, the generation of droplets is fully synchronised with the external perturbation, resulting in droplets of uniform size. This is defined as the synchronised breakup mode (S mode). In the fourth response mode, significant radial broadening can be observed on the jet interface, leading to a pancake-like interface profile. This mode occurs under violent perturbation of the jet flow rate, where fluid bulks rapidly collide in the axial direction, causing the liquid to flatten and extend along the radial direction. According to the experimental study by Meier, Kliipper & Grabitz (Reference Meier, Kliipper and Grabitz1992), this mode is defined as kinematic gathering (K mode). It is notable that the K mode typically exhibits three-dimensional characteristics as the pancake structure evolves until breakup, which lies beyond the scope of our axisymmetric simulation. Therefore, this study only demonstrates the parameter space of the K mode without further analysis of its breakup dynamics.
Figure 3. (a) Four typical response modes of jet breakup under external actuation. The letter symbols I, M, S and K stand for the irregular breakup mode, the multiple breakup mode, the synchronised breakup mode and the kinematic gathering mode, respectively. (b) Phase diagram of different response modes as the frequency
$f$
and the amplitude
$A$
of perturbation vary. The dash-dotted and dashed lines correspond to the theoretical predictions of the synchronised frequency range based on the fully viscous flow (FVF) model and the viscous potential flow (VPF) model of the liquid jet.
Figure 4. (a) Interface profiles of liquid jet as the perturbation frequency
$f$
changes under a fixed perturbation amplitude
$A=2 \times 10^{-4}$
. Temporal evolution of jet tip position at (b)
$f=0$
, (c)
$f=0.04$
, (d)
$f=0.11$
and (e)
$f=0.18$
. The pinch-off positions are also indicated.
Figure 3(b) shows the phase diagram of different response modes. Clearly, the K mode mainly occurs when the perturbation amplitude
$A$
is relatively high (e.g.
$A \sim O(1)$
). Meanwhile, the jet breakup presents the I mode under very low
$A$
(e.g.
$A \sim O(10^{-8})$
), as the external perturbation is too weak to modulate the jet breakup. As the values of
$A$
are in the range
$O(10^{-7} {-} 10^{-1})$
, the response modes are closely related to the perturbation frequency
$f$
. Specifically, there exists a frequency range where the S mode occurs. Below this frequency range, the response mode shifts to the M mode; while beyond this range, the jet breakup presents the I mode, indicating that the external perturbations with high frequencies is unable to dominate the jet breakup. The S mode allows for the generation of uniform droplets with adjustable size and productivity. It can also be observed that as
$A$
increases, the frequency range of the S mode widens significantly. Figure 3(b) also presents the theoretical predictions for the frequency range of the S mode as the amplitude
$A$
changes. In theoretical analyses, both the fully viscous flow (FVF) model and the viscous potential flow (VPF) model of the liquid jet have been used to study the temporal growth rate of perturbation. The details will be provided in § 4.
3.2. Effects of frequency and amplitude
Figure 4(a) shows the interface profiles of the liquid jet as
$f$
changes under a fixed value of
$A=2 \times 10^{-4}$
. For the free liquid jet without external actuation (i.e.
$f = 0$
), the breakup presents the I mode. Under low values of
$f$
(e.g.
$f = 0.03$
and 0.04), the breakup of the liquid jet presents the M mode, with multiple droplets generated in one period of
$1/f$
. It should be noted that the multiple droplets of different sizes will eventually merge together as they evolve downstream (see the case
$f = 0.04$
, for example). The droplet merging dynamics, accompanied by the interaction between droplets and ambient gas in the M mode, will be analysed in detail in § 5. As
$f$
continues to increase (e.g.
$f = 0.06$
, 0.11 and 0.15), the jet breakup presents the S mode, indicating that the size and generation frequency of uniform droplets can be modulated by
$f$
. When
$f$
exceeds the critical value (e.g.
$f = 0.18$
), the jet breakup reverts to the I mode. To illustrate the dynamic characteristics of liquid jet breakup more clearly, figure 4(b–e) presents the temporal evolution of the jet tip position (denoted by
$z_{tip}$
, as sketched in figure 4
a). As the interface breaks up into droplets, the tip position undergoes a sudden jump and the new tip position indicates the jet breakup length. For the I mode at
$f = 0$
and 0.18, as shown in figures 4(b) and 4(e), the pinch-off positions of the interface are random and the jet breakup length varies within a certain range. In the M mode at
$f = 0.04$
(see figure 4
c), multiple pinch-off positions exist within one period of
$1/f$
. In the S mode at
$f = 0.11$
, as demonstrated in figure 4(d), the pinch-off positions remain constant, indicating that the jet breakup is periodic.
Figures 5(a) and 5(b) further present the variations of droplet diameters
$D$
and jet breakup length
$L_j$
under different values of
$f$
. As the surface wave prevents the droplets from forming a standard spherical shape, the values of
$D$
are obtained by calculating the volume of each droplet. The values of
$L_j$
are obtained according to the
$z_{tip}{-}t$
map, as shown in figure 4(b–e). For the unactuated liquid jet (
$f=0$
), which presents the I mode, both the droplet diameter and jet breakup length fall within a certain range, and their average values are also given in the figures (see the symbol accompanied by the error bar). By counting the number of generated droplets in a long time sequence, the average breakup frequency (also referred to as the natural breakup frequency
$f_n$
) can be obtained, with a specific value of 0.107, as indicated in figure 5(a). Under relatively low frequency where the jet breakup presents the M mode, the jet breakup behaviours circulate with multiple droplets generated in one period; therefore, both the values of
$D$
and
$L_j$
are distributed at some discrete values. A continuous increase in
$f$
leads to the jet breakup in the S mode, which results in unique values of droplet diameter and jet breakup length. It should be noted that the natural frequency
$f_n$
invariably falls within the frequency range of the S mode. Since the generation of droplets is fully synchronised with external actuation, the droplet size can be modulated by
$f$
, with the quantitative relationship
$D=(6/f)^{1/3}$
according to the conservation of liquid flow rate. This relationship is confirmed by the numerical results in figure 5(a). Moreover, the jet breakup length is also affected significantly by
$f$
and the minimum length occurs at
$f=0.11$
, which is very close to the natural frequency of jet breakup. The physical reason lies in that the perturbation growth rate of liquid jet reaches its maximum close to the natural breakup frequency, thus causing the earliest breakup of liquid jet as the initial perturbation amplitude of the jet remains unchanged. The theoretical prediction of jet length as
$f$
varies is also given in figure 5(b), as shown by the dashed line. Details of the theoretical analysis will be given in § 4.4.
Figure 5. (a) Droplet diameters at different
$f$
under a constant amplitude
$A = 2 \times 10^{-4}$
, where the natural breakup frequency is denoted by
$f_n$
. (b) Breakup lengths of liquid jets as
$f$
changes at
$A = 2 \times 10^{-4}$
, with the theoretical prediction of jet length also indicated by the dashed line.
Figure 6. (a) Interface profiles of jet as
$A$
changes under a fixed frequency of
$f$
= 0.11. (b) Jet breakup length
$L_j$
with varying
$A$
under a constant
$f$
= 0.11.
We also investigate the influence of actuation amplitude
$A$
on the breakup of liquid jets. It has been clearly observed in figure 3 that external actuation is too weak to dominate the jet breakup under very low
$A$
, leading to the occurrence of the I mode; whereas, under relatively large
$A$
, the external actuation is strong enough to cause the kinematic gathering of the liquid jet, thus resulting in the K mode. As the values of
$A$
are in the range
$O(10^{-7} {-} 10^{-1})$
, the jet breakup can present either M, S or I modes, depending on the value of actuation frequency
$f$
. In this section, we mainly focus on the variation of actuation amplitude
$A$
at a fixed frequency
$f=0.11$
, which is close to the natural breakup frequency. Figure 6(a) shows the jet interface profiles under different values of
$A$
. The jet presents I mode under
$A=0$
and
$2 \times 10^{-8}$
, where the droplet size appears to be non-uniform. At
$A=2 \times 10^{-6}, 2 \times 10^{-4}, 2 \times 10^{-2}$
and
$2 \times 10^{-1}$
, the jet exhibits the S mode, showing periodic behaviour with uniform droplets formed and a consistent breakup length. Only the parameter spaces of the I and S modes are considered here, as the K mode, which occurs under larger values of
$A$
, may present three-dimensional jet breakup characteristics. It is also notable that due to the irregular breakup manner of the I mode, the jet perturbation wavelength at a certain instant could differ from that of the S mode. However, the statistical average value of the perturbation wavelength under the I mode approximates the perturbation wavelength of the actuated jet with
$f=0.11$
. This point can be seen clearly in the statistical droplet size values given in figure 5(a), where the average droplet size of an unactuated jet is approximately equal to that of the actuated jet with
$f=0.11$
. Figure 6(b) further provides the relationship between the jet breakup length
$L_j$
and the actuation amplitude
$A$
. It is notable that the jet breakup length falls within certain ranges with its average value in the I mode, whereas it remains constant in the S mode. The jet breakup length appears to follow the relationship
$L_j \sim$
-ln
$A$
(see the dashed line). The theoretical analysis will also be given in § 4.4.
3.3. Effect of jet velocity
Since the variation of jet velocity can affect the natural breakup characteristics of the liquid jet, it is intended to further influence the response behaviours of the jet to external perturbations. Under the non-dimensional parametric system described in § 2, the jet velocity can be regulated through changing the value of
$\textit{Re}$
at constant values of
$\textit{Oh}$
,
$r_\rho$
and
$r_\mu$
, and this approach has been used by Moallemi et al. (Reference Moallemi, Li and Mehravaran2016). However, this commonly used approach of velocity variation would introduce significant complications for the subsequent quantitative analysis. Specifically, manipulating the jet velocity through the variation of
$\textit{Re}$
would modify both the characteristic velocity (
$\bar U_1$
) and the corresponding characteristic frequency (
$\bar U_1/R$
). According to the definition of
$\textit{Re}$
, the characteristic frequency can be expressed as
$\mu _1\textit{Re}/\rho _1R^2$
. As the characteristic frequency changes with varying jet velocity, a rescale process must be further introduced when analysing the quantitative data related to frequency (e.g. natural breakup frequency, synchronised frequency range), making the theoretical analysis excessively cumbersome. To overcome this limitation, in the current study, the jet velocity is regulated by adjusting the dimensionless mean velocity at the jet inlet while keeping the value of
$\textit{Re}$
constant at the reference state. For example, the dimensionless velocity at the jet inlet was set as
$U_1(t)=1+A \sin (2\pi ft)$
(see (2.6)) in previous sections, whereas in this section, we set the jet velocity as
\begin{eqnarray} U_1(t)=V[1+A \sin (2\pi ft)], \end{eqnarray}
where
$V$
represents the dimensionless mean velocity of the jet, and the product of
$V$
and
$A$
denotes the magnitude of the perturbed velocity. When the jet average velocity equals to unity (i.e.
$V=1$
), (3.1) recovers to (2.6). In this approach, through varying the dimensionless inlet velocity
$V$
without changing the value of
$\textit{Re}$
(= 367), the jet velocity can be adjusted under constant characteristic velocity and frequency, facilitating the further quantitative analyses.
Figure 7. (a) Interface profiles of liquid jets as the dimensionless average velocity
$V$
in (3.1) changes under fixed frequency of
$f$
= 0.11 and amplitude of
$A=2 \times 10^{-4}$
. (b) Breakup lengths of liquid jets as
$A$
changes at
$A = 2 \times 10^{-4}$
and
$f=0.11$
, with the theoretical prediction also indicated by the dashed line. (c) Frequency range of the S mode under different values of
$V$
, where
$f_n$
,
$f_l$
and
$f_u$
denote the natural breakup frequency, the lower critical frequency and the upper critical frequency of the S mode region, respectively. The symbols and the lines denote the numerical results and the theoretical predictions, respectively.
Figure 7(a) shows the interface profiles of liquid jets under different values of
$V$
, where the perturbation amplitude and frequency are fixed at
$A=2 \times 10^{-4}$
and
$f=0.11$
, respectively. The axial coordinates are also given in each graph, presenting the overall increasing tendency of jet breakup length with the continuous increase of
$V$
. It can be observed that when
$V=1$
, 1.5 and 2, the jet presents the S mode, where the jet interface breaks up periodically and the droplets with uniform size and separating distance are generated. The size and separating distance of the uniform droplets increase obviously with the increase of
$V$
. As
$V$
gradually decreases (e.g.
$V=0.5$
), the liquid jet breakup shifts to the I mode, producing droplets with non-uniform size. The result indicates that the external perturbation is unable to modulate the jet breakup with relatively low velocity. While at relatively high jet velocity (e.g.
$V=2.5$
), the jet breakup is found to shift to the M mode, where two droplets are generated in a period of
$1/f$
. The variation of jet length under different values of
$V$
is given in figure 7(b), where the jet breakup length locates with some ranges under the I mode, remains constant under the S mode, and distributes at discrete points under the M mode. The theoretical prediction is also given by the dashed line (described in detail in § 4.4).
To present the jet breakup characteristics in a wide range of jet velocities, figure 7(c) displays the natural breakup frequency
$f_n$
and the frequency range of the S mode (with upper and lower critical frequencies denoted by
$f_u$
and
$f_l$
, respectively) under different jet velocities
$V$
. The symbols represent the numerical results and the lines represent theoretical predictions. It is notable that the statistic error bars for
$f_n$
,
$f_u$
and
$f_l$
in numerical simulations are also given accompanying the average value. The theoretical analysis will be provided later in (4.3) through the linear instability analysis. Clearly, with the increase of jet velocity, both the natural breakup frequency
$f_n$
, and the upper and lower critical frequencies (
$f_u$
and
$f_l$
) of the S mode increase monotonously. The synchronised frequency range is also found to widen notably.
3.4. Development of interface perturbation
The spatial development of interface perturbation is studied quantitatively in this section, as the generation of droplets is directly related to the development of interface perturbation. Without loss of generality, we only consider the liquid jets with fixed axial velocity of unity (2.6) but varying perturbation frequency
$f$
and amplitude
$A$
. It should be noted that, to resolve the interface perturbation more precisely (especially at the upstream position of the jet, where the amplitude of interface perturbation appears to be very small), a much finer mesh size of
$\Delta x=0.02$
is used in the calculation. Figure 8 shows the jet interface profiles accompanied by the spectrums of interface perturbation waves under different response modes. The spectrums of perturbation wavelength are obtained through a statistical approach which is proposed by Zhao, Sprittles & Lockerby (Reference Zhao, Sprittles and Lockerby2019), and this method has been successfully employed in evaluating the dominant instability modes of liquid films (Zhao et al. Reference Zhao, Zhang and Si2023, Reference Zhao, Qiao, Mu, Si and Luo2024). For the liquid jets considered in this study, a discrete Fourier transform is applied to the interface position
$r(z, t)$
at different time instants to obtain the power spectral density of the perturbations
$R(\lambda , t)$
,which is expressed as
\begin{equation} R(\lambda , t)=\sum _{n=0}^{N-1}\, r(x_{n}, t) \, e^{-i \frac {2\pi }{N \lambda } n}, \end{equation}
where
$\{ r(x_{n}) \}_{n=0}^{n=N}$
is the numerical array of the interface profiles. Since
$R(\lambda , t)$
is a complex number, its square root (modulus) can be further evaluated. In this way, the ensemble average value of its square root (denoted by
$R_{rms}$
) can be obtained. It is notable that for each case, at least two hundred interface profiles at different time instants
$t$
are applied for the spectrum analysis. The peak of the spectrum corresponds to the dominant wavelength of interface perturbation.
Figure 8. Interface profile (left column) and wavelength spectrum (right column) for (a) free liquid jet without actuation, (b) S mode under
$A=2 \times 10^{-3}, f=0.15$
, (c) I mode under
$A=2 \times 10^{-2}, f=0.18$
and (d) M mode under
$A=2 \times 10^{-3}, f=0.04$
.
For the free liquid jet without actuation, as shown in figure 8(a), the interface perturbation develops from a very small initial value, which is almost invisible when
$z \leqslant 120$
. As the jet evolves downstream, the interface perturbation grows rapidly, eventually leading to the breakup of the liquid jet. Through spectrum analysis of the jet interface over a time sequence, it is observed that the peak of the spectrum occurs at the wavelength of
$\lambda =8.93$
, which is very close to the wavelength of natural jet breakup (denoted by
$\lambda _n$
, calculated through
$1/f_n$
= 9.35, see figure 5
a) obtained by counting the number of droplets in numerical simulations. Moreover, the spectrum of perturbation lies within a certain range, indicating that the jet breaks up randomly into droplets of non-uniform size (i.e. I mode). For the I mode which occurs under a very small amplitude
$A$
, the interface evolution and perturbation spectrum are similar to the situation of the free liquid jet. It should be noted that as the actuated liquid jet evolves axially with an average velocity of unity, a specific external perturbation with frequency
$f$
corresponds to a certain wavelength of
$1/f$
. For the S mode of jet breakup (see figure 8
b, where
$f=0.15$
and
$A=2 \times 10^{-3}$
), the interface perturbation grows from a relatively small value, with its wavelength almost invariably equal to that of the external actuation
$1/f$
, leading to the generation of droplets with uniform size. In this situation, the spectrum of the interface perturbation is rather narrow and the peak value is almost equal to the corresponding wavelength of the actuation frequency (i.e.
$1/f$
), indicating that external actuation leads to the synchronised breakup of the liquid jet. Figure 8(c) shows the I mode as
$f$
exceeds the critical frequency of the S mode (
$f=0.18$
and
$A=2 \times 10^{-2}$
), where the interface evolution and perturbation spectrum present different characteristics. Specifically, the initial perturbation wavelength on the jet interface is almost equal to
$1/f$
, with its amplitude decreasing as the jet evolves downstream. Near the breakup region, a random perturbation wave with a wavelength approaching
$1/f_n$
can be observed, resulting in the breakup of the liquid jet. Therefore, two peaks can be observed in the spectrum: one corresponding to the wavelength of actuation frequency (
$1/f$
) and the other to the wavelength of natural frequency (
$1/f_n$
). As the perturbation spectrum near the natural frequency lies within a certain range, the droplets formed after the jet breakup show a non-uniform size distribution. For the M mode, which occurs when
$f$
is lower than the frequency range of the S mode, as shown in figure 8(d) (where
$f=0.04$
and
$A=2 \times 10^{-3}$
), the initial perturbation wavelength on the jet interface is almost equal to
$1/f$
. However, the perturbation wave develops into two derivative waves with shorter wavelengths downstream and, eventually, the jet breaks up into two droplets in one period of
$1/f$
. Therefore, three peaks can be observed in the perturbation spectrum. The peak with the largest wavelength approaches the actuation wavelength of
$1/f$
and the other two peaks correspond to the perturbation waves that result in the generation of multiple droplets. The wavelengths of these two derivative waves lie on both sides of the natural wavelength
$1/f_n$
, which is caused by the synergetic influence of both the applied perturbation wave and the random perturbation of the liquid jet under the M mode.
4. Theoretical analysis
4.1. Linear instability analysis
To analyse the mechanisms of jet breakup under external actuation, we perform a linear instability analysis that considers the temporal development of jet perturbations. As shown in figure 9, a liquid jet with axial velocity
$\bar U_1$
emerging into static gas surroundings is considered as the physical model. Similar to the numerical simulations, the dynamic viscosity and density are defined as
$\mu _i$
and
$\rho _i$
, where
$i=1, 2$
represents the liquid jet and the surrounding gas, respectively, and the surface tension is denoted by
$\sigma$
. The axisymmetric cylindrical coordinate system
$(z, r)$
is used to describe the physical model, where
$z$
and
$r$
represent the axial and radial directions, respectively. The characteristic scales in the theoretical analysis are consistent with the numerical simulation (i.e. the jet radius
$R$
and average velocity
$\bar U_1$
are chosen as the characteristic length and velocity, respectively). In this way, the Reynolds number
$\textit{Re}$
, the Ohnesorge number
$\textit{Oh}$
, the density ratio
$r_\rho$
and the viscosity ratio
$r_\mu$
all maintain the same values as in the numerical simulation, with
$\textit{Re}=367$
,
$\textit{Oh}=0.0105$
,
$r_\rho =0.001$
and
$r_\mu =0.025$
, respectively. Similar to the numerical simulation, when the variation of jet velocity is considered in theoretical analysis, the dimensinless jet velocity is adjusted without changing the value of
$\textit{Re}$
.
Figure 9. Theoretical model of a free liquid jet injecting into static gas surroundings in cylindrical coordinates of
$(z, r)$
.
The motions of the liquid jet and the ambient gas are governed by the dimensionless Navier–Stokes equations,
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {u}_i = 0, \\[-28pt] \nonumber \end{align}
\begin{align} \rho _i \left(\frac {\partial \boldsymbol {u}_i}{\partial t}+\boldsymbol {u}_i \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol {u}_i \right) = -\boldsymbol{\nabla }p_i +\frac {1}{\textit{Re}_i}{\nabla} ^2 \boldsymbol {u}_i, \\[9pt] \nonumber \end{align}
where
$\boldsymbol {u}_i$
,
$\rho _i$
and
$p_i$
represent the velocity, density and pressure, respectively, where the subscripts
$i$
= 1 and 2 refer to the liquid jet and the surrounding gas. The Reynolds numbers are defined as
$\textit{Re}_1=\textit{Re}$
and
$\textit{Re}_2=\textit{Re}\boldsymbol{\cdot }r_\rho /r_\mu$
. It is notable that in the numerical simulations, the surface tension force is treated as a body force within the Navier–Stokes equation (see (2.4)). In contrast, in the linear instability analysis, the influence of surface tension is incorporated through the boundary condition at the interface. The theoretical model is solved using the normal mode method, where the flow velocity, pressure and interface position in the Navier–Stokes equations are decomposed into basic quantities and small axisymmetric perturbations with the Fourier form
$\sim e^{\mathrm{i}(k z-\omega t)}$
, where
$k$
is the dimensionless wavenumber, which is inversely proportional to the perturbation wavelength
$\lambda$
. The dimensionless complex frequency
$\omega =\omega _r + \mathrm{i} \omega _i$
consists of
$\omega _i$
, representing the growth rate of perturbations, and
$-\omega _r$
, representing the angular frequency of perturbations. Perturbations with
$\omega _i\gt 0$
are temporally unstable, while those with
$\omega _i \leqslant 0$
are invariably stable. The quantitative relationship between the angular frequency
$-\omega _r$
and the perturbation frequency
$f$
corresponds to
$-\omega _r=2\pi f$
. In this study, the angular frequency
$-\omega _r$
will be used to predict the natural and critical frequencies of liquid jet breakup (see (4.2) and (4.3), respectively), and the growth rate
$\omega _i$
will be used to analyse the jet breakup length, as shown in (4.4). The governing equations and boundary conditions for the perturbations are similar to those of Si et al. (Reference Si, Li, Yin and Yin2009), and the Chebyshev spectral collocation method is employed to solve the eigenvalue problem numerically. Details are provided in Appendix A.
In our study, we also apply the viscous potential flow (VPF) theory (Joseph, Wang & Funada Reference Joseph, Wang and Funada2007) to analyse the growth of perturbations. Compared with the theoretical analysis in Appendix A (referred to as fully viscous flow, FVF), the VPF theory ignores the viscous effects inside the liquid jet and the surrounding gas, and only accounts for the fluid viscosities in the dynamic boundary conditions at the interface. This approach allows for the derivation of an explicit analytical dispersion relation for perturbation growth, with the specific form (Funada, Joseph & Yamashita Reference Funada, Joseph and Yamashita2004)
\begin{eqnarray} && \left [ \frac {I_{0}(k)}{I_{1}(k)} + r_\rho \frac {K_{0}(k)}{K_{1}(k)} \right ] \omega ^2 + 2 \left [ - k \frac {I_{0}(k)}{I_{1}(k)} + i \frac { k^2}{\textit{Re}} \left ( \frac {I_{0}(k)}{I_{1}(k)} - \frac {1}{k} + r_\mu \left ( \frac {K_{0}(k)}{K_{1}(k)} + \frac {1}{k} \right ) \right ) \right ] \omega \nonumber \\ && + \left [ k^2 \frac {I_{0}(k)}{I_{1}(k)} - i \frac {2 k^3}{\textit{Re}} \left ( \frac {I_{0}(k)}{I_{1}(k)} - \frac {1}{k} \right ) -\frac {1}{\textit{We}}k( k^2 - 1) \right ] = 0, \end{eqnarray}
where
$I_0$
and
$I_1$
are the zero-order and first-order Bessel functions, and
$K_0$
and
$K_1$
are the zero-order and first-order modified Bessel functions, respectively. It has been shown that both FVF and VPF theories converge to the inviscid case when the Reynolds number is sufficiently high (e.g.
$\textit{Re} \sim O(10^3)$
or higher). Moreover, the VPF theory retains most of the characteristics predicted by FVF theory for Reynolds numbers ranging from
$O(10^1)$
to
$O(10^2)$
(Funada, Joseph & Daniel Reference Funada, Joseph and Daniel2002; Joseph et al. Reference Joseph, Wang and Funada2007), suggesting that the VPF theory is well suited for the conditions in this study, where
$\textit{Re}=367$
.
4.2. Prediction of the natural frequency
Figure 10. Perturbation growth rate
$\omega _i$
and angular frequency
$-\omega _r$
versus wavenumber
$k$
under
$\textit{Re}=367, \textit{We}=14.9, r_\rho =0.001$
and
$r_\mu =0.02$
, where the fully viscous flow (FVF) model under varying axial jet velocities of
$V=0.5$
, 1, 2 and the viscous potential flow (VPF) model under constant axial jet velocity of unity are considered. The angular frequency (
$-\omega _{\textit{rn}}$
) corresponding to the wavenumber of maximum growth rate (denoted by
$k_m$
) indicates the theoretical prediction of natural breakup frequency
$f_n$
, where
$f_n=-\omega _{\textit{rn}}/(2 \pi )$
.
Figure 10 shows the perturbation growth rate
$\omega _i$
and the angular frequency
$-\omega _r$
of the liquid jet as wavenumber
$k$
varies, where both the FVF model and the VPF model are considered at fixed values of
$\textit{Re}=367, \textit{Oh}=0.0105, r_\rho =0.001$
and
$r_\mu =0.02$
. For the FVF model, the curves under different values of jet mean velocity
$V$
(=0.5, 1 and 2) are also given for comparison. For each curve, it is clear that there exists a cutoff wavenumber beyond which the perturbation growth rate becomes negative, indicating a smallest wavelength of jet instability. Comparing the curves of the FVF model and the VPF model with axial velocity of unity (i.e.
$V=1$
), it is observed that the growth rate of the VPF model remains slightly higher than that of the FVF model, suggesting that the VPF model predicts a more unstable liquid jet. This difference occurs because the VPF model ignores viscosities in the fluid bulk and therefore overestimates the growth rate. A smaller value of
$\textit{Re}$
could lead to a greater discrepancy between the two models (Funada et al. Reference Funada, Joseph and Daniel2002). Comparing the curves with different
$V$
under the FVF model, it is found that the variation of
$V$
only changes the growth rate
$\omega _i$
slightly (see the inset), but mainly affects the perturbation angular frequency
$-\omega _r$
. Specifically, the angular frequency is found to be equal to the product of the wavenumber and the axial jet velocity invariably, indicating that the perturbation propagates with a dimensionless phase velocity of
$-\omega _r/k=V$
. Therefore, if one stands on the axial local framework along with the liquid jet, the interface disturbance only grows temporally and does not propagate upstream or downstream of the jet, similar to the standard temporal instability analysis of liquid jets (Lin Reference Lin2003; Eggers & Villermaux Reference Eggers and Villermaux2008).
According to the linear instability analysis (see figure 10), the wavenumber (
$k_m$
) corresponding to the maximal value of the growth rate curve (
$\omega _{i, \textit{max}}$
) grows fastest among all the perturbation waves, dominating the breakup of the liquid jet. The angular frequency corresponding to the wavenumber
$k_m$
is denoted by
$-\omega _{\textit{rn}}$
, with the quantitative relationship
$-\omega _{\textit{rn}}=k_m \boldsymbol{\cdot }V$
. Therefore, the natural breakup frequency
$f_n$
of the liquid jet can be represented as
$f_n=k_m \boldsymbol{\cdot }V/(2 \pi )$
. For the liquid jet with mean velocity of unity, the FVF model and VPF model give the theoretical results of
$k_m=0.69$
and 0.7, corresponding to the theoretical values of
$f_n=0.109$
and 0.111, respectively. The theoretical predictions obtained by both models agree well with the numerical results shown in figures 5 and 8(a). As the liquid jet velocity gradually changes, the natural breakup frequency can also be easily predicted according to the value of
$k_m$
. The theoretical results under different values of jet velocity have been given in figure 7(c), which present a good agreement with the numerical results.
4.3. Prediction of the synchronised breakup region
Apart from the natural breakup frequency, the frequency range of the synchronised breakup region (i.e. S mode) under different perturbation amplitudes
$A$
can also been given theoretically through comparing the breakup time between an unactuated jet and the forced liquid jet. For a liquid jet with a dimensionless radius of unity, the perturbation on the interface develops with a dimensionless growth rate
$\omega _i$
and an initial magnitude
$\eta _0$
. Under the framework of linear instability, they follow the relationship
\begin{equation} \eta _0 \mathrm{e}^{\omega _i t_b}=1, \end{equation}
where
$t_b$
is the dimensionless jet breakup time, which can be expressed as
\begin{equation} t_b=\frac {1}{\omega _i} \mathrm{ln}\frac {1}{\eta _0}. \end{equation}
Previous research by Moallemi et al. (Reference Moallemi, Li and Mehravaran2016) has indicated the relationship between the initial surface perturbation
$\eta _0$
and the velocity perturbation
$A$
through an energy conservation analysis, with the form (see (14) of Moallemi et al. (Reference Moallemi, Li and Mehravaran2016))
\begin{equation} \dfrac{3}{2}A=\eta _0+(\eta _0-A)\frac {1}{\textit{We}_j}, \end{equation}
where
$\textit{We}_j$
(=
$\textit{We} \boldsymbol{\cdot }V^2$
) denotes the equivalent Weber number of the jet. The equation can be rewritten as
\begin{equation} \eta _0=\frac {1.5\textit{We}_j+1}{\textit{We}_j+1}\boldsymbol{\cdot }A. \end{equation}
Therefore, the relationship between
$t_b$
and
$A$
can be expressed as
\begin{equation} t_b=\frac {1}{\omega _i} \mathrm{ln}\bigg[\frac {\textit{We}_j+1}{(1.5\textit{We}_j+1)\boldsymbol{\cdot }A}\bigg]. \end{equation}
As the perturbation growth rate
$\omega _i$
has been provided in figure 10, the jet breakup time
$t_b$
for different values of
$k$
can be determined. Figure 11 shows the
$t_b{-}k$
curves accompanying with the angular frequency
$-\omega _r$
for the FVF and VPF models under constant jet velocity of
$V=1$
(corresponding to
$\textit{We}_j=\textit{We}=14.9$
), considering the variations of actuation amplitude
$A$
from
$2 \times 10^{-8}$
to
$2 \times 10^{-3}$
. For these cases, the angular frequency
$-\omega _r$
remains equal to that of the wavenumber
$k$
invariably. It should be noted that for the free jet without external actuation (
$A=0$
), the initial perturbation amplitude
$\eta _0$
is too small to be detected, making it impossible to calculate the jet breakup time directly using (4.5). Nevertheless, the dimensionless breakup time of the free liquid jet can be obtained from the jet breakup length using the estimation
$t_b = L_j / V$
. In the numerical simulations, the average value of
$L_j$
is 183 (see figure 6
b), indicating a theoretical value of
$t_b=183$
for an unactuated liquid jet, as shown by the dash-dot-dotted transverse line in figure 11. By comparing the
$t_b{-}k$
curves for external actuation under different
$A$
with the breakup time of the unactuated jet, the wavenumber region where actuation affects the jet breakup can be identified. It is observed that the values of
$t_b$
are invariably larger than 183 when
$A=2 \times 10^{-8}$
, indicating that such a small actuation amplitude has little effect on modulating jet breakup. This result is consistent with the mode diagram in figure 3(b), where jet breakup remains in the I mode under
$A=2 \times 10^{-8}$
. As the perturbation amplitude
$A$
gradually increases, the wavenumber region where
$t_b\lt 183$
corresponds to the perturbation waves that cause the actuated jet to break up earlier than the unactuated jet. Thus, within this wavenumber region, jet breakup can be synchronised with the actuation, resulting in droplet generation under the S mode. By comparing the curves under different amplitude
$A$
, we can determine the wavenumber region for the S mode. The corresponding angular frequency region can be further obtained through the
$-\omega _r \,{\rm versus }\,k$
curve, as sketched in figure 11, where
$-\omega _{\textit{ru}}$
and
$-\omega _{\textit{rl}}$
denote the upper and the lower angular frequencies of the S mode, respectively. The corresponding frequency of
$f_u$
and
$f_l$
can be calculated by
$f_u=-\omega _{\textit{ru}}/(2 \pi )$
and
$f_l=-\omega _{\textit{rl}}/(2 \pi )$
, respectively. In the phase diagram of jet response modes (see figure 3
b), the predicted synchronised frequency range under different perturbation amplitude
$A$
has been given by the lines. Overall, both the FVF model and the VPF model can give a relatively accurate prediction of the synchronised frequency region, with the FVF model appearing to be closer to the numerical results. However, comparing with the FVF model where the growth rate curve is solved through spectral method, the VPF model provides the specific expression for perturbation growth rate and wavelength, as shown by (4.3). The analytical equation allows us to derive the growth rate curve directly, which would facilitate the further calculations of the synchronised frequency range. The simplified form of the VPF model holds important practical significance in engineering applications.
Figure 11. Jet breakup time
$t_b$
and angular frequency
$-\omega _r$
versus wavenumber
$k$
under different perturbation amplitude
$A$
, where the dimensionless jet velocity equals to unity. The solid lines and the dashed lines represent the growth rate curves obtained by the FVF model and VPF model, respectively, while the dash-dot-dotted transverse line in the
$t_b{-}k$
map represents the breakup time of a free liquid jet without external perturbation. The frequencies
$-\omega _{\textit{ru}}$
and
$-\omega _{\textit{rl}}$
indicate the theoretical prediction of upper and lower frequencies of the S mode (
$f_u$
and
$f_l$
), where
$f_u=-\omega _{\textit{ru}}/(2 \pi )$
and
$f_l=-\omega _{\textit{rl}}/(2 \pi )$
, respectively.
It should be noted that the lower boundary of the synchronised frequency region predicted by theoretical analysis shows a significant deviation from the numerical results at relatively high values of
$A$
(e.g.
$A \geqslant 2 \times 10^{-3}$
), as shown in figure 3(b). A qualitative explanation is provided here. On one hand, the growth rate of perturbation is relatively small under low actuation frequencies, making jet breakup more susceptible to random perturbations. As a result, a narrower synchronised frequency range is obtained in numerical simulations compared with theoretical predictions. On the other hand, as
$A$
increases, the jet breakup length decreases, causing the nonlinear stage of perturbation development to play a more significant role in the overall breakup process. Since our theoretical analysis focuses only on the linear instability of jet breakup, greater deviations between theoretical analysis and numerical simulations are expected at relatively high values of
$A$
.
Figure 12. Jet breakup time
$t_b$
and frequency
$-\omega _r$
versus wavenumber
$k$
obtained by FVF model, where the jet velocity
$V$
varies under constant amplitude
$A=2 \times 10^{-4}$
. The dash-dot-dotted transverse lines in the
$t_b{-}k$
map represent the breakup time of unactuated liquid jets with different
$V$
, and the upper and lower frequencies of the S mode are indicated by
$-\omega _{\textit{ru}}$
and
$-\omega _{\textit{rl}}$
, with
$f_u=-\omega _{\textit{ru}}/(2 \pi )$
and
$f_l=-\omega _{\textit{rl}}/(2 \pi )$
, respectively.
The frequency range of the S mode as the jet velocity changes can also been predicted by comparing the breakup time of the unactuated jet and the forced liquid jet. For simplicity, we consider the liquid jet with varying mean velocity
$V$
, but constant actuation amplitude of
$A=2 \times 10^{-4}$
, and the results are shown in figure 12. It is notable that the variation of
$V$
results in the change of
$\textit{We}_j$
(e.g.
$\textit{We}_j$
= 3.7, 14.9 and 59.5 under
$V$
= 0.5, 1 and 2, respectively). As the growth rate of perturbation with different
$V$
only changes slightly (see figure 10), the difference of the breakup time
$t_b$
under different wavenumber
$k$
is mainly caused by the variation of
$\textit{We}_j$
, as shown clearly by (4.8). Similar to the situation considered in figure 11, the breakup time for the unactuated jet under different velocity
$V$
can be obtained through dividing the average jet breakup length
$L_j$
by
$V$
, and the results are given by the transverse lines in the
$t_b{-}k$
map of figure 12. Comparing the results among
$V=0.5$
, 1 and 2, it is observed that the variation of jet velocity only affects the breakup time slightly. The reason can be explained as follows: since the jet propagates with relatively small
$\textit{Oh}$
(=0.0105), the liquid viscosity only plays a tiny role, and the breakup time is decided by the capillary time scale (
${\sim} \sqrt {\rho _1R^3/\sigma }$
). Therefore, the breakup time remains constant theoretically as
$R$
does not change. According to figure 12, the corresponding critical wavenumbers and angular frequencies of the synchronised breakup region for each value of
$V$
can be obtained by the junctions of curves. Through the relationships
$f_u=-\omega _{\textit{ru}}/(2 \pi )$
and
$f_l=-\omega _{\textit{rl}}/(2 \pi )$
, the upper and lower frequencies of the S mode as the jet velocity
$V$
changes can be further predicted. The theoretical predictions have been given by the lines in figure 7(c), showing a good agreement with the numerical results.
4.4. Prediction of jet breakup length
The jet breakup length
$L_j$
can be calculated by multiplying the jet velocity (non-dimensionalised as
$V$
) and the breakup time
$t_b$
, i.e.
\begin{equation} L_j=\frac {V}{\omega _i} \mathrm{ln}\frac {1}{\eta _0}=\frac {V}{\omega _i} \mathrm{ln} \left[\frac {\textit{We}_j+1}{(1.5\textit{We}_j+1)\boldsymbol{\cdot }A}\right]\!. \end{equation}
As
$\textit{We}_j=\textit{We} \boldsymbol{\cdot }V^2$
, (4.9) indicates that the jet breakup length is inversely proportional to the growth rate of perturbation and in positive correlation with the jet velocity
$V$
. Once the jet velocity remains constant (i.e. fixed
$V$
and
$\textit{We}_j$
), the jet breakup length follows the relationship of
$L_j \sim -\mathrm{ln} A$
at fixed
$f$
as the growth rate
$\omega _i$
is decided by
$f$
. Here, we further provide a quantitative comparison between the theoretical predictions and the numerical results as
$f$
,
$A$
and
$V$
change individually. Since the FVF model is more reasonable to the real flow situation and offers a more accurate prediction of perturbation development compared with the VPF model, only the results of the FVF model are used here. As
$f$
changes under constant
$A$
(see figure 5
b, where
$A=2 \times 10^{-4}$
), the corresponding perturbation growth rate
$\omega _i$
can be obtained for different values of angular frequency
$-\omega _r$
and wavenumber
$k$
according to the growth rate curve in figure 10. Therefore, the theoretical predictions of
$L_j$
as
$f$
varies (
$\omega _i$
changes under constant
$V=1$
and
$\textit{We}_j=14.9$
) can be given by the dashed lines in figure 5(b). It is clear that the theoretical results agree well with the numerical simulations near the natural frequency, but diverge significantly as
$f$
deviates from the natural frequency. The reason lies in that as
$f$
diverges from the natural frequency, the perturbation growth rate decreases significantly; thus, the random initial perturbations of the liquid jet have a more pronounced effect on the jet breakup, causing the breakup length to approach that of the free liquid jet, as shown in figure 5(b).
The jet breakup length as
$A$
changes under fixed
$f = 0.11$
is also given in figure 6(b). For the theoretical analysis, a fixed
$f = 0.11$
corresponds to
$-\omega _r=k = 0.69$
, leading to a growth rate of
$\omega _i = 0.089$
, based on the FVF curve in figure 10. Thus, the jet breakup length can be obtained theoretically as
$A$
varies from
$2 \times 10^{-8}$
to
$2 \times 10^{-2}$
under constant
$V=1$
and
$\textit{We}_j=14.9$
, as shown by the dashed line in figure 6(b). Similarly, the jet breakup length under varying values of velocity
$V$
can also be calculated under fixed values of
$f$
and
$A$
, and the theoretical results have been given in figure 7(b). It is clear that the theoretical predictions agree well with the numerical results under the S mode (
$V=1$
, 1.5 and 2), but a certain degree of deviation can be observed under the M mode (
$V=2.5$
and 3), which is similar to the results in figure 5(b). Moreover, as the external actuation fails to modulate the jet breakup under the I mode where
$V=0.5$
, the numerical jet breakup length is obviously larger than the theoretical prediction.
4.5. Further discussions
As the linear instability analysis is able to predict the breakup characteristics (i.e. natural frequency, synchronised frequency region and jet length) of the actuated jet, a further discussion is carried out in this section, aiming to summarise the connection between the theoretical analyses and the numerical simulations, which would also provide guidance for practical applications. In applications, as the medium and velocity of the liquid jet are decided, the key dimensionless parameters (i.e.
$\textit{Re}$
,
$\textit{Oh}$
, density ratio
$r_\rho$
and viscosity ratio
$r_\mu$
) can be calculated precisely. Based on these dimensionless parameters, the linear instability analysis can be carried out, which provides the growth rate/frequency-wavenumber (
$\omega _i/\omega _r{-}k$
) curves of perturbation, where the wavenumber
$k$
is directly related to the actuating frequency of
$f$
. Consequently, the natural frequency of jet breakup can be determined from the peak of this curve, as shown in figure 10. The obtained natural frequency provides guidance for us in applying external disturbances close to this frequency, ensuring the production of uniform droplets in applications. At a given actuating frequency, the jet breakup length can also be further controlled by adjusting the perturbation amplitude, with the corresponding relationship between dimensionless parameters given by (4.9). Additionally, the
$\omega _i{-}k$
curve enables the derivation of the jet breakup time curve (
$t_b{-}k$
) under different actuating amplitudes and jet velocities. By comparing the
$t_b{-}k$
curve with the jet’s free breakup time (which is equal to the division of jet length by velocity), we can determine the synchronised frequency range for different perturbation amplitudes and jet velocities (as shown by figures 11 and 12). This analysis provides guidance for selecting appropriate frequency ranges in practical applications to achieve uniform droplet generation.
5. Droplets merging dynamics under M mode
As multiple droplets generated within one period would eventually merge together under the M mode (see figure 4
a for example), this section analyses the dynamics of droplets merging in detail. Figure 13(a) shows a typical case where two droplets are generated in one period at
$f=0.04$
and
$A=2 \times 10^{-3}$
. The droplets are labelled as ‘D1’ and ‘D2’ based on their sequence of generation, while the merged droplet downstream is labelled as ‘D3’. The temporal evolutions of the mass centre (denoted by
$z_c$
) and the average velocity (denoted by
$V_c$
) of droplets are shown in figures 13(b) and 13(c), respectively. For the convenience of analysis, the zero time instant is set at the moment when the second droplet (i.e. ‘D2’) is generated. It is observed that ‘D2’ gradually approaches and eventually merges with ’D1’, as shown at times
$t_1$
,
$t_2$
and
$t_3$
in figures 13(a) and 13(b). The average velocity of ‘D2’ remains larger than ‘D1’ invariably and the vibrating of droplets caused by surface tension can lead to the pulsation of jet average velocity, as shown in figure 13(c). It is also notable that the average velocities of droplets all remain slightly smaller than the initial jet velocity at the inlet (i.e. 1), which is caused by the non-negligible resistance of the static air surroundings.
Figure 13. (a) Interface profile of droplets in the M mode with
$f=0.04$
and
$A=2 \times 10^{-3}$
, where the two droplets generated within one period are labelled as
$D1$
and
$D2$
, and the merged droplet is labelled as
$D3$
. (b) Temporal evolutions of the droplet centre of mass
$z_c$
. (c) Temporal evolutions of the droplet average velocity
$V_c$
.
Figure 14. (a) Interface profile of droplets in the M mode with
$f=0.03$
and
$A=2 \times 10^{-3}$
, where the three droplets generated within one period are labelled as
$D1$
,
$D2$
and
$D3$
, and the merged droplets are labelled as
$D4$
and
$D5$
. (b) Temporal evolutions of the droplet centre of mass
$z_c$
. (b) Temporal evolutions of the droplet average velocity
$V_c$
.
For the situation where more than two droplets are generated within one period, the dynamics of droplet merging becomes more complex. Figure 14(a) shows a typical case where three droplets are generated in one period under
$f=0.03$
and
$A=2 \times 10^{-3}$
, where the three droplets are labelled as ‘D1’, ‘D2’ and ‘D3’ according to their sequence of generation. The droplet formed by the merging of ‘D1’ and ‘D2’ is labelled as ‘D4’, and the droplet formed by the merging of ‘D4’ and ‘D3’ is labelled as ‘D5’. The temporal evolutions of the mass centre (denoted by
$z_c$
) and the average velocity (denoted by
$V_c$
) of droplets are shown in figures 14(b) and 14(c), respectively. The zero time instant is set at the moment when the last droplet in the period (i.e. ‘D3’) is generated. Combining the results from figures 14(a), 14(b) and 14(c), it is evident that the velocity of ‘D2’ consistently remains higher than that of ‘D1’, leading to the merging of the droplets and the formation of ‘D4’ around the instant of
$t_2$
. Similarly, the velocity of ‘D3’ remains higher than that of ‘D4’, resulting in the merging of the droplets and the final formation of ‘D5’. Similar to the situation shown in figure 13, the average velocities of all droplets remain slightly lower than the initial jet velocity due to air resistance.
Figure 15. Pressure field and streamlines around the droplets under constant
$A=2 \times 10^{-3}$
with (a)
$f=0.04$
and (b)
$f=0.03$
.
The flow fields around the droplets are also analysed to reveal the mechanism of droplet merging. Figure 15(a) shows the pressure field around the droplets under
$f=0.04$
and
$A=2 \times 10^{-3}$
. It is observed that during the downstream propagation of droplets, a high pressure zone occurs near the front edge of droplet ‘D1’, causing the deceleration of ‘D1’. The droplet ‘D2’ with a higher initial velocity gradually approaches ‘D1’, leading to their eventual merging. As for the velocity field close to the droplets, the streamlines are also plotted in figure 15(a). For the convenience of analysis, a translation of coordinates is implemented, where the velocity along the flow direction is subtracted by the average velocity of the liquid jet at the inlet. Therefore, the problem can be shifted to the aerodynamic deformation of liquid droplets in an incoming air flow (from right to left in figure 15
a) with dimensionless velocity of unity. At the leeward of the droplet, a recirculation flow region forms due to flow separation. The low-pressure zone in the recirculation region, combined with the high-pressure zone at the droplet’s front edge, results in the deceleration of droplet. It should be noted that although the streamlines appear to penetrate the droplet interface due to the unsteady flow characteristics, it does not indicate the fluids exchange between the liquid and the gas phases. This is because that the streamlines only depict the directional field of velocity at a frozen instant, without representing the temporal evolution process of fluid elements. Figure 15(b) shows the pressure field accompanying the streamlines (after the same coordinate translation as figure 15
a) under
$f=0.03$
and
$A=2 \times 10^{-3}$
, in which three droplets (denoted as ‘D1’, ‘D2’ and ‘D3’) are generated sequentially within one period. Similarly, a high-pressure zone forms near the front edge of the droplet train, causing the droplets to decelerate as they move downstream. The recirculation flow also occurs at the leeward side of the droplets due to flow separation.
6. Manipulation of droplets merging under M mode
As the M mode results in the coalescence of multiple droplets generated within one period, the manipulation of the droplet merging process is further analysed in this section. The controllable merging of droplets is of practical significance in applications such as the formation of tin droplets with large separation distances in extreme ultraviolet lithography. The large separation distance eliminates the interference between neighbouring droplets during the droplet breakup under strong laser exposure. In this study, we aim to control the dynamics of the M mode by adding a strong pulse (with amplitude
$A_p$
) to the sinusoidal perturbation wave (with dimensionless frequency
$f$
and amplitude
$A$
) every
$n$
(= 2, 3, 4, …) periods, as sketched in figure 16. This additional pulse introduces an equivalent frequency of
$f_p$
(
$=f/n$
). The strategy is designed to first generate a train of uniform droplets via the sinusoidal perturbation wave, then merge
$n$
droplets into one through the addition of the pulse. The merging process is expected to be controlled by the relative amplitude
$A_p/A$
between the additional pulse and the sinusoidal wave.
Figure 16. Sketch of the velocity perturbation at the jet inlet, where a pulse with amplitude
$A_p$
is applied on the sinusoidal perturbation wave (with frequency
$f$
and amplitude
$A$
) every
$n$
periods, resulting in an additional perturbation wave with a frequency
$f_p=f/n$
.
Figure 17. (a) Interface profiles of liquid jets as
$n$
= 2, 4 and 6 under constant
$f=0.16$
,
$A=2 \times 10^{-3}$
and
$A_p=2A$
, corresponding to
$f_p$
= 0.08, 0.04 and 0.027, respectively. In these cases, the growth rate of perturbation with
$f_p$
remains larger than that with
$f$
. (b) Interface profiles of liquid jets under pure sinusoidal perturbation with
$A=4 \times 10^{-3}$
and varying
$f$
.
It should be noted that this strategy to manipulate the dynamics of the M mode is closely related to the growth rates of perturbations corresponding to the specific frequencies of
$f$
and
$f_p$
, respectively. Since the amplitude of the pulse is definitely larger than the sinusoidal perturbation wave, the perturbation growth rate of the pulse frequency should be chosen smaller than that of the sinusoidal wave frequency to ensure the feasibility of this manipulating strategy. Otherwise, the jet breakup is intended to be fully decided by the pulse frequency as the perturbation with pulse frequency
$f_p$
would grow faster than that with frequency
$f$
. To demonstrate the failure of this manipulating strategy, figure 17(a) considers the situation of varying
$n$
(i.e.
$f_p$
) under constant values of
$f=0.16$
,
$A=2 \times 10^{-3}$
and
$A_p=2A$
. For each case, the growth rate of the pulse remains significantly larger than that of the sinusoidal perturbation. For example, the pulses with frequency
$f_p=f/2$
,
$f/4$
and
$f/6$
correspond to the dimensionless growth rate of 0.077, 0.044 and 0.03, respectively, according to the FVF model curve in figure 10, which are much larger than the growth rate of the sinusoidal wave (approaching zero at
$f=0.16$
). It is observed that the jet breaks into uniform droplets for
$n=2$
with an equivalent frequency of
$f_p=0.08$
, and into two and three droplets with equivalent frequencies of
$f_p$
for
$n=4$
and 6, respectively, indicating the failure of this manipulation strategy. In these situations, the jet dynamics is similar to the response to a pure sinusoidal perturbation with frequency of
$f_p$
and amplitude
$A_p$
. For example, under the pure sinusoidal perturbation with amplitude
$A=4 \times 10^{-3}$
(see figure 17
b), the jet breakup exhibits the S mode with uniform droplets generated at
$f=0.08$
, while the jet breakup presents the M mode with two and three droplets generated in one period at
$f=0.04$
and 0.027, respectively. The numerical results are in accordance with the hypothesis discussed previously, as the jet breakup is dominated by the frequency
$f_p$
.
When the growth rate of the pulse is smaller than that of the sinusoidal perturbation wave, the jet breakup dynamics is found to be more complex as both the sinusoidal perturbation and the additional pulse can significantly influence the breakup, depending on the relative amplitude between the pulse and the sinusoidal perturbation. Figure 18(a) shows a group of cases as
$A_p/A$
changes under constant
$f=0.12$
,
$A=2 \times 10^{-3}$
and
$n=4$
, corresponding to a constant
$f_p=0.03$
. According to the growth rate curve (FVF model) of figure 10, the growth rate of the sinusoidal wave with
$f=0.12$
is equal to 0.085, and that of the pulse with
$f_p=0.03$
is equal to 0.034. It is observed that the value of
$A_p/A$
significantly affects the jet breakup dynamics. Under relatively low
$A_p/A$
(e.g.
$A_p/A=2$
), droplet generation is governed by the sinusoidal perturbation wave, as the formation of droplets is synchronised with
$f$
. In this situation, the jet interface profile still resembles that of the S mode, where the disturbance with wavelength
$1/f$
propagates and amplifies downstream. However, the presence of the pulse with
$f_p$
results in a certain degree of a droplet merging as the droplet train evolves downstream. Since the pulse amplitude is relatively low, the final merging of four droplets into one takes a long distance, which is not depicted in the figure. With the increase of
$A_p/A$
(e.g.
$A_p/A=4$
and 8), faster merging of the four droplets within one pulse period is observed, as the pulse plays a more prominent role, demonstrating the feasibility of modulating the merging dynamics of the droplet series. In these cases, the jet interface morphology is supposed to present a distinct superposition of two wavelengths (i.e.
$1/f$
and
$1/f_p$
), thus causing the jet first break into uniform droplets with wavelength
$1/f$
and then coalesce to one large droplet with equivalent wavelength of
$1/f_p$
. However, with a further increase in
$A_p/A$
(e.g.
$A_p/A=10$
and 15), the dynamics of jet breakup becomes different. In these cases, the jet breakup presents the M mode with three droplets generated in one peroid of
$1/f_p$
. The results indicate that the jet breakup is totally dominated by the pulse frequency, as the single sinusoidal perturbation with
$f=0.03$
and the same amplitude of the pulse (i.e.
$A=0.02$
and 0.03) could lead to the M mode with three droplets in one period. The corresponding results are also shown in figure 18(b). In this situation, the wavelength of
$1/f_p$
is supposed to dominate the perturbation of the jet surface over that of
$1/f$
. Overall, we can conclude that the manipulation of the pulse is only effective within a certain range of
$A_p/$
A. A pulse with a very large amplitude results in the domination of the pulse, while a pulse under the critical amplitude can affect the merging distance of the uniform droplets generated by the sinusoidal perturbation.
Figure 18. (a) Interface profiles of liquid jets as
$A_p$
varies under constant
$f=0.12$
,
$A=2 \times 10^{-3}$
and
$n=4$
, corresponding to
$f_p=0.03$
. In these cases, the growth rate of perturbation with
$f_p$
is lower than that with
$f$
. (b) Interface profiles of liquid jets under pure sinusoidal perturbation with
$f=0.03$
and varying
$A$
.
Figure 19. Jet response dynamics with the variations of
$f_p$
(realized by changing the value of
$n$
) and
$A_p/A$
, where the filled circles denote the controllable region where
$n$
droplets merge to one periodically, the hollow diamonds denote the pulse-dominated region.
The number of merged droplets can also be adjusted by changing the value of
$n$
, which alters the equivalent pulse frequency
$f_p$
. For example, with a constant value of
$f=0.12$
, six droplets can be merged together when
$n=6$
, corresponding to
$f_p=0.02$
; five droplets can be merged together when
$n=5$
, which corresponds to
$f_p=0.024$
; and so on. It is found that the critical pulse amplitude for the controllable merging of droplets depends on the value of
$f_p$
, as shown in figure 19. Generally, a lower
$f_p$
(larger
$n$
) leads to a larger critical pulse amplitude for the controllable region. This can be explained by comparing the perturbation growth rates of the sinusoidal wave and the pulse under different values of
$f_p$
. As shown in figure 10 (see the FVF model curve), the growth rate of the sinusoidal wave at
$f=0.12$
equals to 0.085, and the growth rates of the pulse with
$n=2$
, 3, 4, 5 and 6 (corresponding to
$f_p=0.06$
, 0.04, 0.03, 0.024, 0.02) equal to 0.063, 0.044, 0.034, 0.026 and 0.022, respectively. Since the growth rate decreases monotonously with the decrease of
$f_p$
, a larger pulse amplitude is needed to result in the domination of pulse frequency.
7. Conclusions
The forced breakup of liquid jets in static gas surroundings is studied through direct numerical simulations, where the Navier–Stokes equations are solved in combination with a diffuse interface method. A sinusoidal perturbation of jet velocity is applied at the inlet, and the effects of perturbation frequency and amplitude on the evolution of the liquid jet are investigated. Four typical jet breakup response modes are observed, which are the irregular breakup mode, the multiple breakup mode, the synchronised breakup mode and the kinematic gathering mode. A phase diagram of modes transition is provided, showing how they vary with the changes in perturbation frequency and amplitude. Specifically, in the synchronised breakup mode, the breakup of the liquid jet produces a train of droplets with uniform size and separation distance, and the droplet size can be adjusted by varying the frequency. The natural frequency of jet breakup is located within the frequency range of the synchronised breakup region. The jet breakup length reaches its minimum when the external frequency approaches the natural frequency, and increasing the actuation amplitude decreases the jet breakup length monotonically. The liquid jet velocity is found to modulate the jet breakup length and also increase the natural breakup frequency, and the upper and lower critical frequencies of the synchronised region. A spectral analysis of interface perturbation development is conducted, revealing the characteristics of different jet breakup response modes. The linear instability analysis based on the fully viscous model and the viscous potential flow model is also performed, providing theoretical predictions of the natural frequency and the frequency range where uniform droplets can be produced. By combining the linear instability analysis with scaling analysis, the variation tendencies of jet breakup length as a function of perturbation frequency, amplitude and jet velocity can be predicted. The merging dynamics of droplets under the multiple breakup mode is also analysed, showing that for the multiple droplets generated in one period, the upstream droplet has a larger velocity compared with the downstream droplet, which results in the the merging of droplets. The high pressure zone forms at the front of droplet train due to air resistance, which results in the deceleration of droplets. Furthermore, by adding a pulse to the sinusoidal perturbation waveform, the merging dynamics of multiple droplets can be modulated actively by adjusting the equivalent frequency and amplitude of the pulse. In the future work, a combination of a series of frequencies can also be applied to the jet, which is expected to bring in richer jet breakup dynamics due to the competition of perturbation waves with different frequencies and growth rates. In summary, this numerical investigation can provide guidance for the experimental study of the modulation of liquid jet breakup, which is expected to contribute to the on-demand generation of droplets in practical applications.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2025.10833.
Funding
This work was supported by the National Natural Science Foundation of China (Nos. 12388101 and 12272372), Strategic Priority Research Program of the Chinese Academy of Sciences (XDB0910100), Youth Innovation Promotion Association CAS (Nos. 2018491, 2023477) and USTC Research Funds of the Double First-Class Initiative (YD2090002020).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Equations for full viscous flow model
For the full viscous flow model applied in the linear instability analysis, the normal mode method is used, in which all the variables such as velocity
$\boldsymbol {u}_i(r,t)$
, pressure
$p_{i}(r,t)$
and interface position
$r_{j}(r,t)$
in (4.2) are split into a basic quantity and a perturbation part:
\begin{equation} \boldsymbol {u}_i=\boldsymbol {U}_i+\tilde {\boldsymbol {u}}_i, p_i=P_i+\tilde {p}_i, r_j=1+\tilde {\eta }, \end{equation}
in which the perturbation components (
$\tilde {\boldsymbol {u}}_{i}$
,
$\tilde {p}_i$
,
$\tilde {\eta }$
) are expressed in Fourier form as
$e^{\mathrm{i}(k z-\omega t)}$
in cylindrical coordinates
$(r, z)$
, with the form:
\begin{equation} \begin{array}{ll} \tilde {\boldsymbol {u}}_{i}=(\hat {u}_{i}(r), \hat {v}_{i}(r)) e^{\mathrm{i}(k z-\omega t)},\quad \tilde {p}_i=\hat {p}_i(r) e^{\mathrm{i}(k z-\omega t)},\quad \tilde {\eta }=\eta e^{\mathrm{i}(k z-\omega t)}, \end{array} \end{equation}
where
$\hat {u}_{i}(r)$
and
$\hat {v}_{i}(r)$
are the amplitudes of velocity perturbations in the
$z$
and
$r$
directions, respectively. Here,
$\hat {p}_{i}(r)$
is the amplitude of pressure perturbation and
$\eta$
is the perturbation amplitude on the jet interface. Additionally,
$k$
is the dimensionless wavenumber, which is the reciprocal of the axial perturbation wavelength
$\lambda$
,
$\omega =\omega _r+\mathrm{i} \omega _i$
is the dimensionless complex frequency, where
$\omega _i$
represents the growth rate of perturbations.
The perturbations of all quantities are assumed to be very small, which allows us to expand the Navier–Stokes equations and only keep the first-order perturbation terms. Therefore, substituting these perturbations into the governing equations and neglecting the high-order terms, we can get the linearised governing equations:
\begin{align} \frac {\mathrm{d}\hat {v}_i}{\mathrm{d}r}+\frac {\hat {v}_i}{r}+\mathrm{i}k\hat {u}_i=0, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {1}{\textit{Re}_i}\left[\frac {\mathrm{d}^2\hat {v}_i}{\mathrm{d}r^2}+\frac {1}{r}\frac {\mathrm{d}\hat {v}_i}{\mathrm{d}r}-\left(k^{2}+\frac {1}{r^2}+\mathrm{i}k \textit{Re}_{i}U_{i}\right)\hat {v}_i \right]-\left(\frac {1}{r_\rho }\right)^{\delta _{i2}}\frac {\mathrm{d}\hat {p}_i}{\mathrm{d}r}+\mathrm{i}\omega \hat {v}_i=0, \\[9pt] \nonumber \end{align}
\begin{align} \frac {1}{\textit{Re}_{i}}\left[\frac {\mathrm{d}^2\hat {u}_i}{\mathrm{d}r^2}+\frac {1}{r}\frac {\mathrm{d}\hat {u}_i}{\mathrm{d}r}-(k^{2}+\mathrm{i}k \textit{Re}_{i}U_{i})\hat {u}_i \right]-\frac {\mathrm{d}U_{i}}{\mathrm{d}r}\hat {v}_i-\left(\frac {1}{r_\rho }\right)^{\delta _{i2}}\mathrm{i}k\hat {p}_i+\mathrm{i}\omega \hat {u}_i=0. \end{align}
The velocity and the pressure at the symmetry axis must satisfy the consistency conditions:
\begin{equation} \hat {v}_1=\frac {\mathrm{d}\hat {u}_1}{\mathrm{d}r}=\frac {\mathrm{d}\hat {p}_1}{\mathrm{d}r}=0. \end{equation}
At the interface
$r_j=1+\eta e^{\mathrm{i}(k z-\omega t)}$
, the continuity of velocity, the kinematic boundary condition, and the force balance on the tangential and normal directions should be satisfied, i.e.:
\begin{align} \hat {v}_1=\hat {v}_{2}, \hat {u}_1+\frac {\mathrm{d}U_1}{\mathrm{d}r}\eta =\hat {u}_{2}+\frac {\mathrm{d}U_{2}}{\mathrm{d}r}\eta , \\[-28pt] \nonumber \end{align}
\begin{align} -\mathrm{i}\omega \eta =\hat {v}_j-\mathrm{i}k U_1\eta , \\[-28pt] \nonumber \end{align}
\begin{align} \hat {p}_{2}-\frac {2r_\mu }{\textit{Re}}\frac {\mathrm{d}\hat {v}_{2}}{\mathrm{d}r}-\hat {p}_1-\frac {2}{\textit{Re}}\frac {\mathrm{d}\hat {v}_1}{\mathrm{d}r} =\frac {1}{\textit{We}}(1-k^2)\eta , \\[-28pt] \nonumber \end{align}
\begin{align} \mathrm{i}k\hat {v}_1+\frac {\mathrm{d}\hat {u}_1}{\mathrm{d}r}+\frac {\mathrm{d}^2 U_1}{\mathrm{d}r^2}\eta =N \left(\mathrm{i}k\hat {v}_{2}+\frac {\mathrm{d}\hat {u}_{2}}{\mathrm{d}r}+\frac {\mathrm{d}^2 U_{2}}{\mathrm{d}r^2}\eta \right)\!. \\[9pt] \nonumber \end{align}
The boundary conditions at infinity (
$r \rightarrow \infty$
) are
\begin{equation} \hat {v}_2=\frac {\mathrm{d}\hat {u}_2}{\mathrm{d}r}=\frac {\mathrm{d}\hat {p}_2}{\mathrm{d}r}=0. \end{equation}
In the calculation, this condition can be satisfied by employing a finite and sufficiently large distance of the gas phase.
The governing equations and the boundary conditions form an eigenvalue problem, which is solved using the Chebyshev spectral collocation method. A MATLAB code was developed for this purpose and more details on the calculation process can be found from Si et al. (Reference Si, Li, Yin and Yin2009).
Open access
