# Bag Film Breakup of Droplets in Uniform Airflows

K. Tang1<EMAIL_ADDRESS>T.A.A. Adcock1  W. Mostert1  1Department of
Engineering Science, University of Oxford, Oxford OX1 3PJ, UK

###### Abstract

We present novel numerical simulations investigating the bag breakup of liquid
droplets. We solve the two-phase incompressible Navier-Stokes equation using
the adaptive mesh refinement (AMR) technique, with interface reconstruction by
a volume-of-fluid (VOF) method. We first show good agreement with theoretical
predictions and clarify the viscous effect on drop deformation before the
onset of bag breakup. Next, in numerical simulations of the bag breakup
regime, a bag film is developed at late time and is susceptible to spurious
mesh-induced breakup, which has prevented previous studies from reaching grid
convergence of fragment statistics. We therefore adopt a recently developed
manifold death (MD) algorithm which artificially perforates thin films once
they reach a prescribed critical thickness independent of the grid size. We
show grid convergence of fragment statistics when utilising the MD algorithm,
and analyse the fragment behaviour and physical mechanisms contributing to bag
breakup including ligament breakup, node detachment and rim destabilisation.
We also examined the viscous effects on the bag breakup process and fragment
statistics, where a non-monotonic dependency of the average diameter of bag
film fragments on the liquid viscosity is found. These results highlight the
utility of this numerical method for two-phase simulations involving thin-film
breakup, and pave the way for physics-based numerical investigations into
spume generation at the air-sea interface.

##  1 Introduction

Liquid atomisation refers to the process where a bulk volume of liquid
disintegrates into fragments featuring various sizes and shapes (Guildenbecher
et al., 2009; Pairetti et al., 2018). The fragments generated are described as
sprays, which are involved in many natural and industrial processes, including
ocean-atmosphere interactions (Veron, 2015; Erinin et al., 2019),
precipitation and rain-drop dynamics (Villermaux & Bossa, 2009; Jalaal &
Mehravaran, 2012; Veron, 2015), combustion of liquid propellant in aerospace
applications (Young, 1995), pharmaceutical spray generation (Mehta et al.,
2017) and pathogen transmission (Bourouiba, 2021; Kant et al., 2022), etc.
More specifically, it has recently been found that the atomisation of small-
scale sea surface perturbations dominates ocean spume generation under extreme
wind conditions, producing large droplets with typical sizes of $10^{2}\sim
10^{3}\,{\rm\mu m}$ (Troitskaya et al., 2017, 2018). In this size range the
currently-available sea-spray generation functions (SSGF), crucial for
calculations of air-sea momentum and heat exchange in earth-system modelling,
show large range of scatter (Veron, 2015). However, since the physics
governing the fragmentation of bag films have not yet been firmly established,
their influences on SSGFs have been difficult to quantify. Improving this
understanding is the primary motivation of the present work.

Two stages of liquid atomisation have been identified within literature,
namely the primary and secondary atomisation. Sheets, ligaments and droplets
are stripped from a bulk fluid during primary atomisation, which further
decompose until stabilising capillary effects take over during secondary
atomisation (Pairetti et al., 2018). Secondary atomisation is typically
modelled by the droplet aerobreakup problem characterised by the interaction
between an initially spherical droplet with density $\rho_{l}$, viscosity
$\mu_{l}$ and diameter $d_{0}$, and an ambient gas flow with density
$\rho_{a}$, viscosity $\mu_{a}$ and uniform velocity $U_{0}$ (Guildenbecher et
al., 2009). Based on the previously mentioned physical properties, together
with surface tension $\sigma$ at the liquid-gas interface, four non-
dimensional controlling parameters have been proposed using Buckingham’s Pi
Theorem (see e.g. Table 1 in (Jalaal & Mehravaran, 2014)):

$We\equiv\frac{\rho_{a}U_{0}^{2}d_{0}}{\sigma},\quad
Oh\equiv\frac{\mu_{l}}{\sqrt{\rho_{l}d_{0}\sigma}},\quad\rho^{*}\equiv\frac{\rho_{l}}{\rho_{a}},\quad\mu^{*}\equiv\frac{\mu_{l}}{\mu_{a}}.$
(1)

Among these, $We$ and $Oh$ are respectively the Weber and Ohnesorge number
quantifying the ratio of inertial to capillary and viscous to capillary
forces, and $\rho^{*}$ and $\mu^{*}$ are respectively the density and
viscosity ratios of the liquid and gas phase.

Within literature, various droplet aerobreakup regimes have been observed
where the droplet shows different deformation patterns, and the transition
thresholds between these regimes have traditionally been delineated using $We$
and $Oh$ (Yang et al., 2017; Zotova et al., 2019; Marcotte & Zaleski, 2019),
although some recent works have shown that the density ratio $\rho^{*}$ may
also play an important role (Yang et al., 2016; Marcotte & Zaleski, 2019; Jain
et al., 2019). $Oh$ has been reported to influence the transition thresholds
only when exceeding a critical value of 0.1 (Hsiang & Faeth, 1995); and as
$We$ increases, the breakup becomes more violent and vibrational, bag, multi-
mode (bag-stamen), sheet-thinning and catastrophic breakup regimes are
observed in succession (Jalaal & Mehravaran, 2014; Kékesi et al., 2014).
Alternatively, based on the governing hydrodynamic instability involved in the
process, the four breakup regimes mentioned above can be re-grouped into two
major categories: Rayleigh-Taylor piercing (RTP) and shear-induced entrainment
(SIE) (Theofanous, 2011). However, despite the extensive amount of related
work, the underlying physics governing the transient drop deformation in each
regime are still largely unclear. Furthermore, the empirical transition
criteria proposed so far are often contradictory (Theofanous, 2011; Yang et
al., 2017), with the notable exception of a consensus that the critical Weber
number beyond which bag breakup initiates is ${We}_{c}=11\pm 2$ when $Oh<0.1$
(Guildenbecher et al., 2009; Yang et al., 2017).

In the bag breakup regime, the initially spherical droplet first flattens and
forms a disc, whose centre is then blown downstream and inflates into a hollow
bag attached to a toroidal rim. The time it takes for the drop to reach
breakup $\Delta t_{d}$ typically falls within the range of $\tau\leq\Delta
t_{d}\leq 2\tau$, where $\tau\equiv\sqrt{\rho_{l}/\rho_{a}}d_{0}/U_{0}$ is the
characteristic deformation time proposed by Nicholls & Ranger (1969). The
swollen bag first ruptures near the centre, triggering expansion of holes on
the surface of the bag and eventually bursting into a large number of
fragments, which is then followed by the breakup of the remnant toroidal rim
into smaller amounts of fragments (Chou & Faeth, 1998; Guildenbecher et al.,
2009; Opfer et al., 2014). Droplet bag breakup and its associated fragment
size and velocity distribution functions are of specific interest as they bear
a strong resemblance to the previously mentioned bag-mediated fragmentation of
small-scale sea-surface perturbations under extreme wind conditions
(Troitskaya et al., 2017, 2018).

Droplet aerobreakup involves a complex interplay of aerodynamic, capillary and
viscous effects that is still poorly understood (Jain et al., 2015). The
prevalent theoretical understanding is that hydrodynamic instabilities,
particularly Kelvin-Helmholtz (KH) and Rayleigh-Taylor (RT) instability, play
an important role in the aerobreakup process (Guildenbecher et al., 2009;
Theofanous, 2011; Theofanous et al., 2012; Jackiw & Ashgriz, 2021). KH
instability occurs at the interface between two different streams of fluid
with different velocities and densities (Kundu et al., 2012). In the context
of large-$We$ droplet aerobreakup, it governs the SIE breakup category
(Theofanous, 2011), and is typically found near the drop periphery where the
relative velocity between the liquid and gas phases is the largest
(Gorokhovski & Herrmann, 2008; Jalaal & Mehravaran, 2014). However, due to
strong capillary effects, KH instability is unable to influence droplet
deformation in the bag breakup regime (Theofanous et al., 2012; Jalaal &
Mehravaran, 2014). RT instability occurs when a corrugated interface
separating fluids with different densities undergoes constant acceleration
(Zhou et al., 2021), and is hypothesised to cause interfacial perturbation
growth on the windward surface of the droplet. The wave number of such
perturbations determines whether the droplet undergoes oscillatory
deformation, bag breakup or multi-mode breakup (Yang et al., 2017). However,
instability theories have difficulty in accounting for the viscous effects
(Jalaal & Mehravaran, 2014), flow dynamics prior to drop flattening, and
finite thickness and peripheral boundary of the flattened disc (Jackiw &
Ashgriz, 2021). Alternatively, some works highlight the influence of the
internal flow within the droplets on the deformation process (Guildenbecher et
al., 2009; Villermaux & Bossa, 2009; Jackiw & Ashgriz, 2021; Obenauf & Sojka,
2021). The internal flow model compensates for the drawback of the RT
instability model in predicting early-time drop deformation; however, this
approach is somewhat empirical and cannot account for the complex interaction
between wake vortices and drop surface (Marcotte & Zaleski, 2019). The late-
time breakup behaviour, on the other hand, is delineated into a bag-film
rupturing event, and the fragmentation of the remnant rim at a later time. The
bag film rupture occurs more rapidly and produces much smaller fragments
compared with the remnant rim breakup, and is thus more difficult to capture
experimentally (Guildenbecher et al., 2009). It has only recently been
clarified (Jackiw & Ashgriz, 2022) that the major pathways leading to bag
fragmentation are the destabilisation and collision of hole rims as they
recede over the curved bag and experience centripetal acceleration. The size
and velocity distribution data of aerobreakup fragments are still scarce (Zhao
et al., 2011), and it remains unclear how each breakup mechanism contributes
to the fragment statistics given the large span in time and length scales.
Furthermore, the effects of the $Oh$ value on the bag breakup fragment
statistics and the underlying physical mechanisms still remain largely unclear
(Jackiw & Ashgriz, 2022).

The earliest research on droplet aerobreakup are mostly experimental, where
the droplet breakup behaviour is recorded and analysed using shadowgraphs,
high-speed cameras and particle image velocimetry (PIV) (Hsiang & Faeth, 1992;
Guildenbecher et al., 2009; Jalaal & Mehravaran, 2012; Radhakrishna et al.,
2021). Thanks to the recent development of computational power, numerical
studies have provided a way to investigate atomisation phenomena and gain
insight into fundamental mechanisms that are otherwise difficult to achieve
experimentally (Gorokhovski & Herrmann, 2008; Ling et al., 2015). However,
serious challenges are also present for computational studies on droplet
aerobreakup, including reaching numerical convergence at large density ratio
$\rho^{*}$ (Zotova et al., 2019; Marcotte & Zaleski, 2019); high computational
cost of fully resolving small-scale fragmentation processes in two-phase
turbulence simulations at high $We$ values (Gorokhovski & Herrmann, 2008;
Jalaal & Mehravaran, 2014; Shinjo, 2018), where the smallest droplet size may
be much less than the Kolmogorov scale (Shinjo, 2018). There is also potential
need of ensemble averaging when fragments produced from an individual
realisation is not sufficient for obtaining statistically meaningful results
(Mostert et al., 2022). In particular, as the Navier-Stokes equations do not
describe the physical mechanisms that control topological changes at phase
boundaries, thin films are subject to uncontrolled numerical perforation when
their thickness approaches the minimum grid size (Chirco et al., 2022). As a
result, the fragment statistics are dependent on grid sizes (Jackiw & Ashgriz,
2022), and numerical convergence with respect to bag fragment statistics has
not previously been obtained to our knowledge. It is therefore of paramount
importance to improve the grid resolution level and make the onset of breakup
independent of the grid size, even though the exact physical mechanism
initiating the breakup events remains elusive (Kant et al., 2022). A few
attempts have been made to improve the numerical resolution of fragmentation
or coalescence phenomena. Among these, Coyajee & Boersma (2009) first proposed
a modified VOF scheme that utilises multiple marker functions for different
fluid interfaces to minimise spurious coalescence on coarse meshes.
Afterwards, Zhang et al. (2019) built a topology-based numerical scheme which
automatically refines grid cells containing the liquid film bordered by two
adjacent bubble interfaces. Finally, Chirco et al. (2022) developed an
algorithm that randomly perforates thin films once their thickness reduces to
a prescribed critical value independent of the grid size. This algorithm is
controllable via a set of tuning parameters and has been shown to improve grid
convergence behaviour for various two-phase problems including droplet
aerobreakup.

We present results of novel multiphase direct numerical simulation of droplet
bag breakup using both axisymmetric and fully three-dimensional
configurations. We conduct axisymmetric simulations to study pre-breakup
deformation dynamics, and three-dimensional studies to investigate the breakup
dynamics and acquire fragment data. Our study is structured as follows. We
present in §2.1 the configuration of our problem and the parameter space we
explore, and then introduce the numerical method in §2.2. We analyse the
axisymmetric simulation results in §3 and compare them with previous
theoretical predictions, focusing on the early-time deformation period where
Jackiw & Ashgriz (2021) predicted a constant spanwise growth rate (§3.1), and
the film-thinning period immediately before bag breakup where an exponential
decay model of film thickness is available (Villermaux & Bossa, 2009) (§3.2).
We then investigate the breakup of the bag film based on the three-dimensional
simulation results, where we first show grid convergence of fragment
statistics using the MD algorithm (Chirco et al., 2022) (§4.1). We then
analyse the size and velocity distributions, and provide an overview for the
breakup mechanisms leading to bag disintegration in §4.2. Afterwards, we track
and reconstruct the evolution of individual fragments, and study the
dependence of their ejection velocity, lifetime and oscillation patterns in
§4.3. Finally, we investigate the influence of $Oh$ values on the breakup of
bag films (§4.4). We conclude the study in §5 with some remarks on future
work.

##  2 Formulation and methodology

###  2.1 Problem description

The flow configurations for axisymmetric and three-dimensional simulations are
shown in Figs. 1a and 1b, respectively. For both axisymmetric and three-
dimensional simulations, a stationary liquid droplet with diameter $d_{0}$,
density $\rho_{l}$ and viscosity $\mu_{l}$ is placed close to the left
boundary, surrounded by an initially quiescent gas phase with density
$\rho_{a}$ and viscosity $\mu_{a}$. The domain width $D$ is set as $10d_{0}$
and $15d_{0}$ for axisymmetric and three-dimensional simulations,
respectively, so as to eliminate the influence of finite domain size on the
aerobreakup process. A zero-gradient velocity boundary condition is applied at
the right boundary and a uniform incoming velocity $U_{0}$ is imposed on the
left boundary, while no-penetration conditions are applied at the other domain
boundaries. This velocity initialisation results in an impulsive acceleration
of the droplet at the first time step, and induces a flow field satisfying
both the incompressible constraint and the conservation of linear momentum
(Jalaal & Mehravaran, 2014; Marcotte & Zaleski, 2019).

(a)

(b)

Figure 1: Sketches showing the initial configurations of axisymmetric (a) and
three-dimensional (b) droplet aerobreakup simulations. The axis of symmetry is
located at the bottom in (a).

As discussed in §1, the problem is defined by four non-dimensional parameters,
namely the Weber number $We$, the Ohnesorge number $Oh$, the density ratio
$\rho^{*}$ and the viscosity ratio $\mu^{*}$. Since we are interested in air-
water systems, $\rho^{*}$ and $\mu^{*}$ are set as 830 and 55, respectively,
following the earlier work of Pairetti et al. (2018). We vary $We$ between 12
and 25 in our axisymmetric simulations, while in our current 3D simulations we
fix it at 15. In the meantime, $Oh$ is varied between $10^{-4}$ and $0.075$,
which allows for a comprehensive investigation of viscous effects on bag
breakup.

###  2.2 Numerical method

We use the open-source Basilisk numerical library (Popinet, 2019) to solve the
Navier-Stokes equations for two-phase incompressible, immiscible and
isothermal flows, which are written in the following variable-density form,

$\displaystyle\nabla\cdot\bm{u}=0,$ (2)
$\displaystyle\rho\left(\frac{\partial\bm{u}}{\partial
t}+\bm{u}\cdot\nabla\bm{u}\right)=-\nabla
p+\nabla\cdot\left[\mu(\nabla\bm{u}+\nabla\bm{u}^{T})\right]+\sigma\kappa\delta_{s}\bm{n}.$
(3)

Equations (2) and (3) are respectively the continuity and momentum equation,
where $p$ is the fluid pressure. Surface tension effects are incorporated in
the volumetric form $\sigma\kappa\delta_{s}\bm{n}$ within Equation (3), where
$\sigma$ is the surface-tension coefficient, and $\kappa$ and $\bm{n}$ are
respectively the local curvature and normal vector on the interface. The Dirac
delta $\delta_{s}$ is non-zero only on the interface, indicating the local
concentration of surface tension effects (Popinet, 2009, 2018).

The geometric volume-of-fluid (VOF) method is applied in Basilisk to
reconstruct the interface and minimize the parasitic currents induced by
surface tension (Popinet, 2018), which solves the following advective
equation,

$\frac{\partial f}{\partial t}+\bm{u}\cdot\nabla f=0,$ (4)

where $f$ is the VOF function that distinguishes the liquid and gaseous
phases, taking the value of 1 and 0 in the former and latter respectively. For
modelling of surface tension effects, $\delta_{s}\bm{n}$ in Equation (3) is
approximated as the gradient of the VOF function $\nabla f$ using an
adaptation of Brackbill’s method (Brackbill et al., 1992; Popinet, 2009), and
the curvature $\kappa$ is calculated by taking the finite-difference
discretisation of the derivatives of interface height functions (Popinet,
2009). The quad/octree-based AMR scheme based on the estimation of local
discretisation errors of the VOF function gradient $\nabla f$ and flow
velocity $\bm{u}$ is adopted so as to reduce the computational cost at high
resolution levels $L$, which is defined using the minimum grid size,

$\Delta=\frac{D}{2^{L}}.$ (5)

In the bag breakup regime, the onset of fragmentation is preceded by the
inflation of bag structure whose thickness reduces considerably over time.
While the mechanism responsible for the puncture of the bag film has been
extensively discussed (Lhuissier & Villermaux, 2012; Chirco et al., 2022; Kant
et al., 2022), in VOF simulations this is initiated when the local thickness
of the bag decreases to the finest grid size (Ling et al., 2017), causing the
initiation time of breakup and the size of the finest fragments to be grid-
dependent. To circumvent this unphysical and numerically uncontrolled
phenomenon, we adopt the manifold death (MD) algorithm recently developed by
Chirco et al. (2022), which artificially perforates thin films once their
thickness decreases to a prescribed critical value independent of the grid
size. This enables grid convergence to be reached in the fragment size
distributions and related quantities Chirco et al. (2022). This is realised in
the Basilisk framework by first computing quadratic moments of the VOF colour
function $f$ on grid cells with a given signature level $L_{\rm sig}\leq L$.
The signs of the computed quadratic moments indicate the local shape of the
interface. If a film with thickness not larger than $h_{c}\equiv 3D/2^{L_{\rm
sig}}$ is detected, the algorithm randomly creates cubic cavities on the
ligament by directly setting the value of $f$ to that of the other phase with
a probability $p_{\rm perf}$. While the total fluid mass is changed when holes
are created on thin films, the MD algorithm is able to minimise this side
effect by creating cavities with minimum sizes that allow for Taylor-Culick
expansion, and limiting the maximum number of holes perforated at every
iteration. Further discussion, and details of the parameters used for the MD
algorithm in our study are supplied in §4.1.

Before the formation of thin bag films and their subsequent breakup, the
smallest length scale in the aerobreakup problem is the thickness of the
viscous air boundary layer $\delta$ around the droplet, through which momentum
diffuses from the surrounding airflow into the droplet and drives its
deformation. Batchelor’s estimation of the defining length scale of the
droplet $d_{0}$ yields $\delta\sim d_{0}/\sqrt{Re}$, where
$Re\equiv\rho_{a}U_{0}d_{0}/\mu_{a}$ is the freestream Reynolds number. For a
typical droplet in the bag breakup regime, characterised by Weber and
Ohnesorge numbers $We=15$, $Oh=10^{-3}$, this corresponds to $\delta\sim
1.2\times 10^{-2}d_{0}$. The recommended criteria of $\delta/\Delta\geq 2$
(Mostert & Deike, 2020) then sets the grid resolution level $L\geq 12$ for
simulations with domain size $D=15d_{0}$. The highest grid resolution level we
set in our simulations is $L=14$, at which the droplet contour in our
axisymmetric simulations has reached grid independence. The numerical
convergence of fragment statistics will be discussed in detail in §4.1.

Finally, the droplet radius $R_{0}$, incoming flow velocity $U_{0}$, dynamic
flow pressure $p_{0}\equiv\rho_{l}U_{0}^{2}$ and the characteristic
deformation time $\tau$ introduced in §1 provide the natural reference scales
for the length, mass and time quantities that appear in Eqs. (2) and (3), and
will be used to non-dimensionalise the numerical results in the remainder of
this study unless otherwise specified.

##  3 Pre-breakup deformation dynamics

Before the onset of bag breakup, the shape of the deforming droplet remains
largely axisymmetric, although the wake region may have become fully turbulent
and three-dimensional. Many previous numerical aerobreakup studies therefore
conducted axisymmetric simulations for a parametric study. In this section, we
present our axisymmetric results to provide an overview of the pre-breakup
deformation characteristics of the droplet, while also validating our
simulation results by comparing with available analytic models.

###  3.1 Early-time deformation

We first discuss the initiation period of aerobreakup, defined by Jackiw _et
al._ as $0\leq t\leq T_{i}$, where $T_{i}$ is the time when the droplet
reaches its minimal streamwise thickness (Jackiw & Ashgriz, 2021). To provide
an overview of the early-time droplet deformation process characterised by
spanwise flattening, we first present in fig. 2 the droplet contours extracted
from our axisymmetric simulations at various instants within $0\leq t/\tau\leq
0.8$ for two different Ohnesorge numbers $Oh$ of $10^{-3}$ and $10^{-2}$,
while the same Weber number $We=15$ is set for both cases. The radial profile
is shown with $y=0$ as the axis of symmetry.

(a)

(b)

Figure 2: Early-time development of droplet contours for axisymmetric
simulations with Ohnesorge number $Oh=10^{-3}$ (a) and $10^{-2}$ (b), while
the Weber number $We=15$. The axis of symmetry is at $y=0$.

It is found that during the early deformation stage, the windward surface of
the droplet continues moving downstream and pushing liquid to the drop
periphery, leading to the gradual spanwise flattening of the droplet. In the
meantime, a dimple develops on the windward surface that moves downwards and
eventually evolves into a crater on the axis of symmetry for $Oh=10^{-3}$, as
shown in fig. 2a. The leeward side of the droplet remains relatively
stationary after an initial movement to the left. In contrast, fig. 2b shows
that the increase of $Oh$ to $10^{-2}$ postpones the dimple formation on the
windward surface significantly, which only begins to appear at $t/\tau=0.8$.
Previous works have attributed the spanwise flattening of the drop to the
aerodynamic pressure difference between the frontal stagnation point and the
equatorial periphery (Jackiw & Ashgriz, 2021), which drives the internal flow
within the droplet against the restoring effects of surface tension (Marcotte
& Zaleski, 2019). The airflow quickly separates from the leeward surface,
creating a re-circulation region with low pressure which induces little
movement at the leeward interface (Jain et al., 2015). Formation of similar
dimple structures on the windward surface also be observed in fig. 1 of
Marcotte & Zaleski (2019) within the Weber number range of $11.3\leq We\leq
24$ corresponding to bag and bag-stamen breakup.

We briefly examine whether the dimple is caused by RT Instability developing
on the windward surface due to wind acceleration. Li et al. (2019) predict a
critical instantaneous Bond number $Bo_{c}\equiv\rho_{l}\alpha
R_{0}^{2}/\sigma=11.2$, beyond which the windward surface is destabilised.
Here $\alpha$ is the instantaneous acceleration of the liquid droplet. For a
droplet with $We=15,\,Oh=10^{-3}$, our results show $Bo=0.57$ at $t/\tau=0.4$
when the dimple is first observed in fig. 2a, much smaller than the threshold
value of 11.2 predicted by Li et al. (2019). Taking into account that the
liquid is being primarily pushed from the frontal stagnation point to the
upper-left tip around the time of dimple formation ($t/\tau\sim 0.4$ in fig.
2a), together with the $We$ range where it is observed in (Marcotte & Zaleski,
2019), it is likely that the dimple formation is caused by the capillary
pinching effects against fluid influx, and should therefore be viewed as a
precursor of later rim formation.

Recently, Jackiw & Ashgriz (2021) found experimentally that there exists an
early period featuring constant growth rate of the maximum spanwise radius of
the droplet $R$, and proposed the following model for its prediction,

$\dot{R}=\frac{R_{0}T_{bal}}{4\tau^{2}}\left(a^{2}-\frac{128}{We}\right),$ (6)

where $a$ is the axial stretching rate near the frontal stagnation point and
approximated as $a\simeq 6$ in Jackiw & Ashgriz (2021); $\tau$ is the
characteristic deformation time introduced in §1; and $T_{bal}$ is the time
when a constant streamwise deformation rate is reached, taken as $0.125\tau$
according to the experimental results of Jackiw & Ashgriz (2021). We note that
this model is derived assuming ellipsoidal or cylindrical droplet shape and a
balance between aerodynamic and capillary forces during deformation, which
leads to a purely radial internal velocity profile that cancels out the
viscous effects.

(a)

(b)

Figure 3: Measured droplet spanwise growth rate compared with Eq. (6) (Jackiw
& Ashgriz, 2021). Evolution of instantaneous spanwise growth rate $\dot{R}$ at
various $We$ and $Oh=10^{-3}$ (a) and various $Oh$ with $We=15$ (b) are
plotted; and the numerical results are normalised using Eq. (6).

We now further investigate (6) using our axisymmetric numerical results.
Figures 3(a) and (b) show the influence of Weber number $We$ and the Ohnesorge
number $Oh$ on the instantaneous spanwise growth rate $\dot{R}$, respectively.
For the small $Oh$ value of $10^{-3}$, fig. 3a indicates that the spanwise
growth rate $\dot{R}$ reaches a plateau with relatively small variations
around $t=0.3\tau$, where the prediction of (6) matches qualitatively with the
measured $\dot{R}$ values. For cases at $Oh=10^{-3}$, this period of constant
$\dot{R}$ ends around $t=0.55\tau$, after which the spanwise growth rate
$\dot{R}$ reaches a peak aound $t=0.6\tau$ and then decreases, indicating a
deviation from model (6) absent in the analyses of Jackiw & Ashgriz (2021). On
the other hand, fig. 3b suggests that as $Oh$ increases beyond $2.5\times
10^{-3}$, the late-time peaking of $\dot{R}$ gradually attenuates, while the
match with (6) is improved and maintained for longer periods of time, which is
particularly interesting as (6) is derived based on inviscid flow assumptions
and cannot account for viscous influences. Note that Jackiw & Ashgriz (2021)
tested droplets for which $Oh=2.7\times 10^{-3}$; our numerical results are
therefore consistent with their experiment.

Returning to fig. 2a suggests that during the period $0.3\tau\leq t\leq
0.55\tau$ when the constant growth rate $\dot{R}$ is observed, the liquid is
being pushed from the windward surface to the upper-left tip, where the
maximum spanwise radius is reached. However, at $t=0.6\tau$ when the peaking
behaviour is observed, a bulge appears downstream to the upper-left tip and
causes a location shift where the maximum spanwise radius $R$ is reached. This
bulging behaviour is virtually absent when $Oh=10^{-2}$, as shown in fig. 2b,
where the periphery of the droplet contour only flattens over time.

(a)

(b)

Figure 4: Flow fields near the tip of a droplet with $We=20,\,Oh=10^{-3}$ (a)
and $We=15,\,Oh=10^{-2}$ (b) when the peaks in $\dot{R}$ are reached. For both
cases $Oh=10^{-3}$. The non-dimensional times at which (a) and (b) are taken
are respectively $t/\tau=0.62$ and 0.66.

To provide insights into physical mechanisms governing the peak in the
spanwise growth rate observed for low Oh values in fig. 3, we plot in fig. 4
the pressure distribution and streamlines near the drop periphery, when the
peaks in the spanwise growth rate $\dot{R}$ are reached in fig. 3b for
$We=20,\,Oh=10^{-3}$ and $We=15,\,Oh=10^{-2}$. It can be seen that the
surrounding gas flow separates from the droplet surface at the upper-left tip,
creating attached recirculating vortices in its wake with low pressure and
slow fluid motion (Marcotte & Zaleski, 2019; Jain et al., 2019), where the
bulges are located. The pressure difference in the surrounding flow between
the frontal stagnation point and the recirculating region drives the internal
flow within the droplet from the windward surface to the periphery.
Furthermore, the peaks in $\dot{R}$ observed when $Oh\leq 0.005$ are
associated with the formation of a high-pressure region at the bulge on the
droplet periphery, as can be seen in fig. 4a, which is caused by surface
tension and decelerates the flow into the bulge. Further development of the
bulge leads to an increase in the local capillary pressure, which causes the
decrease in $\dot{R}$ after the peak. Notably, the droplet contour in fig. 4b
at $Oh=10^{-2}$ lacks craters at the axis of symmetry and bulges at the
periphery, and therefore more closely resembles the cylindrical geometry
assumed in (Jackiw & Ashgriz, 2021), which may explain why the match with the
inviscid model (6) is improved as $Oh$ is increased.

(a)

(b)

Figure 5: (a): Evolution of air ($p_{a}$) and liquid pressures ($p_{w}$) on
either side of the droplet interface as functions of the interfacial arc
length $l$; (b): axial airflow velocity $u_{z}$ on the axis of symmetry as a
function of the distance to the windward stagnation point of the droplet $z$.
The values of $We$ and $Oh$ are respectively 15 and $10^{-3}$.

We further investigate the distribution patterns of flow pressure and velocity
in the vicinity of the droplet, and their association with prediction (6) of
Jackiw & Ashgriz (2021). Figure 5a shows the air- and liquid-phase pressure on
either side of the drop surface as functions of the arclength $l$ along the
axisymmetric droplet contour. It can be seen that at very early time
$(t/\tau=5\times 10^{-2})$, the air-phase pressure profile closely follows the
sinusoidal potential-flow solution for $l/R_{0}\leq 1.25$, which corresponds
to the windward face of the drop, whereas the profile at $l/R_{0}>1.25$
deviates from the potential-flow solution due to flow separation,
characterised by a second minimum around $l/R_{0}=2.5$. We also note that at
$t/\tau=5\times 10^{-2}$ the shape of the liquid-phase pressure profile bears
strong resemblance to its airphase counterpart, with a nearly-uniform upshift
due to the constant capillary pressure difference $2\sigma/R_{0}$. As the
droplet flattens over time, the air-phase pressure profile on the windward
surface increases and the first minimum moves to the left, deviating from the
potential-flow solution. In the meantime, the change in liquid-phase pressure
for $l/R_{0}\leq 0.75$ is relatively small, and the air- and liquid-phase
pressure profiles cross over each other at $l/R_{0}=0.75$ and $t/\tau=0.4$,
signalling the dimple formation on the windward surface as the local radius of
curvature reaches infinity. It is also noted that the minimum of the liquid-
phase pressure profile around $l/R_{0}=1.75$ observed at $t/\tau=0.4$ becomes
a maximum at $t/\tau=0.6$, which corresponds to the bulge formation observed
in fig. 4a that leads to the deviation from (6) (Jackiw & Ashgriz, 2021).

Figure 5b shows the airphase axial velocity measured on the axis of symmetry
as a function of the distance to the windward stagnation point of the drop
$z/d_{0}$, where the slope of the curves corresponds to the axial stretching
rate $a$ used in model (6). It is first observed that the axial velocity value
at $z/d_{0}=0$ increases gradually over time, which is because the measuring
point is located in the airphase boundary layer attached to the accelerating
droplet. The axial stretching rate $a$ is found to gradually decrease from 6
and approach $\pi/4$, which are the extreme values corresponding to spherical
and pancake drop shapes as noted in Jackiw & Ashgriz (2021). The decrease of
$a$ corresponds to the air-phase pressure increase on the windward surface as
observed in fig. 5a via the following equation:

$p_{g}(r)-p_{g}(0)=-\rho_{g}\frac{a^{2}U_{0}^{2}}{8d_{0}^{2}}r^{2}.$ (7)

Consequently, we conclude that our numerical results reproduce the prediction
(6) of Jackiw & Ashgriz (2021) that there exists a period characterised by
constant spanwise radius growth rate $\dot{R}$. Furthermore, we find that the
later deviation from (6) is characterised by a peak in $\dot{R}$, which is
caused by the capillary deceleration of liquid influx into the drop periphery
that causes the formation of a bulge at low $Oh$ values. The increase in $Oh$
eliminates the bulge and the frontal crater on the droplet surface, causing
the droplet to acquire a nearly cylindrical shape which is one of the
underlying assumptions by Jackiw & Ashgriz (2021) when deriving (6), hence the
better match with their model.

###  3.2 Film drainage and onset of bag breakup

When the droplet deforms into a disc at the end of the early-time deformation
period, corrugations develop on their frontal surfaces which have generally
been considered as RT perturbation waves (Yang et al., 2017). This appears in
fig. 6a at $t/\tau=0.95,\,We=15$, in the form of waves on the windward surface
of the droplet. Later on, thick rims are observed to form at the drop
periphery due to capillary pinching effects, which extract liquid from the
drop centre and contribute to the formation of bag films near the axis of
symmetry. Subject to the aerodynamic pressure difference between their frontal
and leeward surfaces, these films further bulge out from the rim (Jackiw &
Ashgriz, 2021) and cause exponential growth of streamwise bag length before
breakup.

Figure 6a shows non-uniform profiles of the bag thickness $h$, featuring a
neck where a local minimum in $h$ is reached and the film breakup eventually
occurs. The neck moves outwards radially at $We=15$, leaving a thickening
remnant stamen structure developing at the axis of symmetry (Marcotte &
Zaleski, 2019). Figure 6b, on the other hand, shows the deformed drop contours
at $t/\tau=1.73$ for various $Oh$ values; and it can be seen that as $Oh$
increases, the neck becomes less obvious as the distribution of bag film
thickness becomes more uniform.

(a)

(b)

Figure 6: (a): Evolution of axisymmetric droplet contours with $We=15$ and
$Oh=10^{-3}$, where the axis of symmetry is at $y=0$. (b): droplet contours at
$t/\tau=1.73$ with various $Oh$ values and $We=15$.

It has been argued that the breakup of bag films is due to an RT instability
peculiar to thin films rather than a finite-time singularity of the Navier-
Stokes equations (Villermaux & Bossa, 2009). Assuming inviscid flow and
uniform bag thickness, Villermaux & Bossa (2009) derived the following
exponential decay model for the film thickness $h$,

$h(t)\sim d_{0}e^{-\lambda t},$ (8)

where the exponential decay rate $\lambda$ is given as $4/\tau$. We compare
our numerical results with the predictions of model (8) in fig. 7. The neck
thickness $h$ is calculated by measuring the minimum distance between the
windward and leeward surfaces of the deformed drop contour over a time period
of $t_{b}-0.87\tau\leq t\leq t_{b}$, where $t_{b}$ is the time when film
breakup is detected. Logarithmic scale is used for the $y$ axis to facilitate
comparison of the decay rate $\lambda$.

It can be seen that for the $We$ and $Oh$ range presented in fig. 7, the
exponential decay rate $\lambda$ is initially a constant close to the
prediction of (8). This phase, which features a constant thickness decay rate,
roughly corresponds to the period of rim development prior to the ‘bulging’ of
bag films, which is shown in fig. 6. However, the decay rate increases as the
film continues thinning and approaches breakup, similar to the result of Kant
et al. (2022), which becomes more significant as $We$ and $Oh$ decrease. Most
notably, at $We=12$ the thinning rate $a$ continuously increases close to the
onset of breakup, which suggests that an exponential decay law in the form of
(8) does not fully capture the underlying physics for film drainage with
strong surface tension.

(a)

(b)

Figure 7: The development of neck thickness $h$ over time, measured from
simulations with various $We$ with $Oh=10^{-3}$ (a) and various $Oh$ with
$We=15$ (b). Villermaux and Bossa’s prediction (8) is plotted for comparison.

Developing new theoretical models in place of (8) whose predictions match
better with the late-time neck drainage behaviour observed within the $We$ and
$Oh$ range of interest is out of the scope of the current work. However, we
note briefly that the drainage behaviour of bag films under the influence of
aerodynamic pressure difference observed here bears resemblance to the
drainage of liquid films between a free air-water surface and a buoyancy-
driven air bubble (Pigeonneau & Sellier, 2011; Kočárková et al., 2013; Guémas
et al., 2015), where film drainage models are developed based on lubrication
assumptions (see Section 4.2 of (Magnaudet & Mercier, 2020) and references
therein for more detailed discussions). More specifically, Pigeonneau &
Sellier (2011) also showed a deviation from exponential decay of bubble film
thickness under asymptotically large surface tension, which is ascribed to a
finite-time singularity and contrasts with the thin-film RT instability
mechanism proposed by Villermaux & Bossa (2009). However, the major difference
between the bag and the bubble film drainage problem lies in the location of
the neck. For the drainage of bag films, the neck can be formed some distance
away from the axis of symmetry due to a competition of inertia between the
outer rim and the inner stamen (Marcotte & Zaleski, 2019), which complicates
theoretical modelling due to additional difficulties in predicting time-
varying neck locations. In contrast, bubble film drainage always occurs on the
axis of symmetry, as the thin bubble film is connected to an infinitely large
and quiescent liquid domain.

##  4 Breakup of Bag Films

In this section, we analyse our three-dimensional simulation results during
the period when the bag film undergoes disintegration to form small fragments,
and the axisymmetric flow assumption completely breaks down. It is noted that
the later stage when the surrounding main rim collides with the receding
remnant bag is not covered, which will be reserved for future work (Jackiw &
Ashgriz, 2022).

###  4.1 Grid convergence for fragment statistics

While the physical mechanisms responsible for the onset of liquid film breakup
in general remain an active research topic, with various candidates proposed
including chemical or thermal inhomogeneities (Kant et al., 2022), Marangoni
effects (Lhuissier & Villermaux, 2012) or presence of surface contamination
(Néel & Villermaux, 2018), it has been argued that for bag films under normal
acceleration, thickness modulations arise across the film due to the RT
instability, resulting in perforation when the perturbation amplitude becomes
comparable to the film thickness $h$ (Villermaux & Bossa, 2009; Jackiw &
Ashgriz, 2021).

Numerically, the perforation and subsequent fragmentation of thin films in
droplet breakup problems has historically been challenging to represent in a
physically consistent manner. Given that perforation involves a topological
change in the air-water interface, numerical studies have usually employed
interface-capturing techniques such as the geometric VOF approach employed
here (Tang et al., 2021; Jain et al., 2015, 2019); such methods can represent
interfacial topological change without a need for extensive special treatment.
However, such techniques also tend to suffer from an unphysical and
numerically uncontrolled perforation and fragmentation mode in thin films,
which moreover compromises numerical convergence of the statistics of the
resulting fragment populations (Chirco et al., 2022). In this phenomenon, when
the film thickness approaches the local mesh size, it begins to destabilise,
generating a large range of small fragments without a well-defined
fragmentation mechanism; fig. 8a shows a qualitative illustration of this
phenomenon in our own simulations. This phenomenon also appears in images of
droplet breakup in Jain et al. (2015); unfortunately, information on numerical
convergence of fragment statistics is not supplied in that study.

The manifold death (MD) algorithm constructed and implemented by Chirco et al.
(2022) into Basilisk aims to bypass this spurious mode of fragmentation. This
algorithm detects and artificially perforates thin films periodically by
removing liquid mass once their thickness decreases to a prescribed critical
value. The key point is that VOF breakup is circumvented because the film is
perforated at a thickness greater than what is required for VOF breakup to
occur. A limitation of the MD method is that it removes mass from the droplet
in order to generate the hole, which disturbs the momentum and mass
conservation properties of the VOF scheme as implemented in Basilisk. We find
in practice that sufficiently few holes are formed in these simulations that
this mass loss is insignificant (see below).

The holes that result from the MD perforation mechanism resemble those
appearing in experiments (such as Lhuissier & Villermaux (2013); Jackiw &
Ashgriz (2022); Kant et al. (2022)) - see fig. 8b. The algorithm also affords
considerable user control over the frequency and location, for example, of the
perforations. In the present study, we control perforation frequency using a
probabilistic approach which is scaled to be independent of parallelisation,
resolution, and the calling interval of the algorithm, in order to minimise
the number of free parameters governing the perforation problem. The
appropriate choice of probability to match the efficiency of hole generation
in experimental studies, such as those seen in Lhuissier & Villermaux (2012);
Vledouts et al. (2016), is a complicated problem which is left for future
work. Our aim in this study is instead to establish numerical convergence and
verify various aspects of the resulting fragmentation process with experiment
and theory, which in our knowledge has not been established in previous
numerical studies.

All grid convergence test cases we conduct in this section are reloaded from a
single three-dimensional simulation snapshot with an intact bag at
$We=15,\,Oh=10^{-3}$, and the perforation probability and calling interval for
the MD algorithm are set as $p_{\rm perf}=5.7\times 10^{-5}$ and $\Delta
t_{c}=0.5d_{0}/U_{0}$, with grid level $L=12,\,13,\,14$ and signature level
$L_{\rm sig}=11,\,12,\,13$ (see §2 for their definition). The fragment
statistics are collected and output at fixed time intervals until the bags
have fully disintegrated, and then post-processed to obtain time- and
ensemble-averaged data. It is noted that as the test cases run without using
the MD algorithm are deterministic in the sense that VOF breakup appears
repeatably in the same locations of the thin film, while producing a large
number of very small fragments, only one realisation is completed at each $L$
value; whereas the MD algorithm introduces randomness in the perforation
location and subsequent fragment formation, while also reducing the number of
fragments, and therefore multiple realisations are completed for a given set
of $L$ and $L_{\rm sig}$ to generate sufficient total number of fragments for
ensemble-averaging. A full list for the configurations of the three-
dimensional numerical simulations is available in Table 1.

$We$ | $Oh$ | $L$ | $L_{\rm sig}$ | Realisation No. | Category
---|---|---|---|---|---
15 | $10^{-3}$ | 12 | N/A | 1 | Convergence - VOF
15 | $10^{-3}$ | 13 | N/A | 1 | Convergence - VOF
15 | $10^{-3}$ | 14 | N/A | 1 | Convergence - VOF
15 | $10^{-3}$ | 13 | 11 | 10 | Convergence - MD
15 | $10^{-3}$ | 13 | 12 | 5 | Convergence - MD
15 | $10^{-3}$ | 13 | 13 | 3 | Convergence - MD
15 | $10^{-3}$ | 14 | 13 | 7 | Convergence - MD
15 | $10^{-3}$ | 14 | 13 | 3 |  $Oh$ Study
15 | $10^{-4}$ | 14 | 13 | 2 |  $Oh$ Study
15 | $5\times 10^{-4}$ | 14 | 13 | 1 |  $Oh$ Study
15 | $10^{-3}$ | 14 | 13 | 3 |  $Oh$ Study
15 | $5\times 10^{-3}$ | 14 | 13 | 3 |  $Oh$ Study
15 | $10^{-2}$ | 14 | 13 | 3 |  $Oh$ Study
15 | $2.5\times 10^{-2}$ | 14 | 13 | 3 |  $Oh$ Study
15 | $5\times 10^{-2}$ | 14 | 13 | 3 |  $Oh$ Study
15 | $7.5\times 10^{-2}$ | 14 | 13 | 1 |  $Oh$ Study
Table 1: List of ensemble realisations for three-dimensional numerical
simulations carried out in this work, where the drop Weber and Ohnesorge
numbers $We$ and $Oh$, the grid and signature levels $L$ and $L_{\rm sig}$,
the number of individual realisations, and the purpose for using the ensemble
data (the grid convergence study for §4.1 and §4.2, or the $Oh$ effect study
in §4.4) are indicated.

(a)

(b)

(c)

Figure 8: Effect of the MD algorithm on the bag breakup behaviour at grid
level $L=12$ and 13 for $We=15,\,Oh=10^{-3}$. (a)-(c): Simulation snapshots
showing fragmenting bag films at $t/\tau=1.909$ without (a) and with
artificial perforation (b,c). The grid resolution level is $L=12$ for (a,b)
and 13 for (c), while the MD signature level for (b,c) is $L_{\rm sig}=12$.

We first demonstrate the effects of the MD algorithm on the breakup behaviour
of the bag film in fig. 8. The film rupture behaviour is qualitatively
different with and without application of the MD algorithm. Figure 8a is a
snapshot for a simulation case run without using the MD algorithm, featuring
numerous small-scale irregular corrugations and ligament breaking on the bag,
which reflect the uncontrolled nature of VOF breakup. Figures. 8b and 8c show
that the MD algorithm is able to create holes on the bag film in a controlled
manner, and reduce the influence of VOF breakup on the bag dynamics. These
holes created by the MD algorithm feature well-defined bordering rims that
recede over the bag and create surface capillary waves ahead of them, and may
collide with one another and form long stretching liquid bridges (Agbaglah,
2021), which are not observed with VOF breakup. Note that fig. 8b still shows
some VOF breakup behaviour which is absent in fig. 8c. This reflects the fact
that, even though the film is perforated when the film thickness reaches the
order of $3\Delta_{\rm sig}$, the perforation probability and the rate of hole
expansion are sufficiently low such that there are regions of the film that
continue to thin down to the order of $\Delta$, where VOF breakup begins. In
fig. 8b, $\Delta=\Delta_{\rm sig}$, so that VOF-induced fragmentation still
appears, but this can be further minimized by choosing $L>L_{\rm sig}$, such
as in fig. 8c.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9: Time- and ensemble-averaged size (left column) and speed (right
column) probability distribution functions of aerobreakup fragments obtained
from simulations without using the MD algorithm (upper row), from an
individual realisation (middle row) and from ensemble-averaged data across
various realisations with the MD algorithm applied (lower row) at various grid
resolution and signature levels. The log-normal fit of the experimental data
of Guildenbecher et al. (2017) is included in (e) for comparison. For all test
cases, $We=15$ and $Oh=10^{-3}$.

Figure 9 further compares the grid convergence behaviour for the size
distribution of bag fragments without and with the MD algorithm applied.
Figure 9a shows the fragment size distribution functions obtained from
individual simulations without application of the MD algorithm. It can be seen
that while the distributions have similar shapes at various grid levels $L$,
i.e. featuring large number densities of small fragments near the minimum grid
size, followed by a fall-off at large fragment sizes, there is no clear
indication of the distribution functions reaching grid convergence. In
particular, it is observed that as $L$ increases, the entire size distribution
shifts to smaller sizes. In contrast, fig. 9c presents the fragment size
distribution functions obtained from individual realisations within the range
of $13\leq L\leq 14$ and $11\leq L_{\rm sig}\leq 13$ when the MD algorithm is
used. While more scatters in the size distribution functions are seen when
compared with fig. 9a due to smaller amounts of fragments produced, we no
longer observe the shift to small fragment sizes for the distribution tail
with $d\geq 8\Delta_{\rm sig}$, which is the range for well-resolved fragments
as observed in (Chirco et al., 2022). Figure 9e further presents the ensemble-
averaged size distribution functions obtained with the MD algorithm applied,
which features much smaller range of scatter and shows clearly that grid
convergence is reached for statistics of fragments satisfying $d\geq
8\Delta_{\rm sig}$. Thus together figs. 9a, 9c and 9e establish that the lack
of grid convergence in the no-MD case (fig. 9a) is attributable not to the
scatter of individual realisations, but specifically to numerically
uncontrolled VOF-breakup. We also include the log-normal model of
Guildenbecher et al. (2017) in fig. 9e, which is fitted in that study to the
size distribution of bag-breakup fragments with $We=13.8$. An excellent
agreement between this log-normal model and the tail of our size distribution
is observed, which further validates our numerical results within the size
range of $d\geq 8\Delta_{\rm sig}$. While we leave for future work the
detailed investigation of possible differences in fragmentation mechanisms
between our present results and the experiments of Guildenbecher et al. (2017)
and Jackiw & Ashgriz (2021, 2022), the present remarkable agreement with
experimental data at larger fragment sizes suggests that the upper tails of
the size distribution do not depend on whatever these differences may be. We
note also that Guildenbecher et al. (2017) express most confidence in this
upper tail of the distribution.

It is also noted that in fig. 9e, the fragment statistics are not fully
converged for $d\leq 8\Delta_{\rm sig}$. This is probably because the
fragments within this range are primarily formed due to the breakup of liquid
ligaments, especially the smallest fragments near the grid size, which are
most likely the satellite drops produced from the capillary breakup of
corrugated slender ligaments (Pal et al., 2021). These are controlled by the
grid level $L$ rather than the signature level $L_{\rm sig}$, as the geometry-
specific MD algorithm only targets thin liquid films in 3D simulations and do
not act in the stretch-induced breakup of liquid ligaments. Our numerical
results for $d\leq 8\Delta_{\rm sig}$ also deviate from the log-normal
function fit of Guildenbecher et al. (2017), which may be due to multiple
factors including the difference in $We$, the presence of additional flow
perturbations in experiments, and possibly resolution limits in experimental
equipment, as exemplified by a comparison performed in fig. 9 of Guildenbecher
et al. (2017). Finally, we remark that the influence of MD on mass
conservation is minimal as the loss of liquid mass incurred by the MD
algorithm does not exceed 0.023% for $t/\tau\leq 2.18$ at $L=13$ and 14.

Figures 9b, 9d and 9f also show the average speed $\bar{v}$ of fragments as
functions of the fragment diameter $d$ obtained from simulations without and
with the MD algorithm applied. Similar to the size distribution functions, the
shapes of the distribution of $\bar{v}$ clearly indicate grid convergence for
the tail constituted by the well-resolved fragments with $d\geq 8\Delta_{\rm
sig}$ in fig. 9f, which is not observed in fig. 9b where large scatters across
various grid levels are present. Interestingly, in fig. 9f, fragments with
diameter $d\geq 8\Delta_{\rm sig}$ show little variation in the average speed,
which appears to be a constant independent of the values of $L$ and $L_{\rm
sig}$; whereas a peak can be observed within the range of $\Delta_{\rm
sig}\leq d\leq 8\Delta_{\rm sig}$ which is grid-dependent, along with the
large increase of the speed of tiny fragments close to $\Delta_{\rm sig}$
which approaches the freestream velocity $U_{0}$. While it is clear therefore
that the production mechanisms of droplets satisfying $d\leq 8\Delta_{\rm
sig}$ may not be grid-converged, the resulting dynamics of these small
droplets turn out to be well-resolved; more detailed analysis establishing
this will follow in §4.2.

In summary, our results in this section demonstrate that the application of
the MD algorithm helps to establish grid convergence of fragment size and
speed statistics for well-resolved fragments with diameter $d\geq 8\Delta_{\rm
sig}$, which is not achieved when VOF breakup is dominant. Based on these
results, all following three-dimensional studies of bag film breakup are
conducted at $L=14$ and $L_{\rm sig}=13$.

###  4.2 Mechanisms leading to bag fragmentation

In this section, we further analyse the fragment statistics obtained from our
grid convergence tests run at $L=14$ and $L_{\rm sig}=13$, to provide insight
into the shapes of the size and distribution functions observed in §4.1, and
the physical mechanisms governing the formation of fragments and their
subsequent evolution patterns. These choices of $L$ and $L_{\rm sig}$,
together with the MD parameters specified in §4.1, allow the creation of only
a few holes on the bag film, which are not only enough to avoid the onset of
VOF breakup, but also preserves abundant film breakup phenomena including rim
recession, collision and destabilisation behaviour that would otherwise be
hard to recover with more holes created, where rim collision would dominate.
It is noted that at $L_{\rm sig}=13$, thin films with thickness around
$3D/2^{L_{\rm sig}}=1.1\times 10^{-2}R_{0}$ are susceptible to MD nucleation,
which is one order of magnitude larger than the experimental value of
$h/R_{0}=1.4\times 10^{-3}$ as found by Jackiw & Ashgriz (2022). Thicker
liquid films lead to smaller Taylor-Culick velocity values, which reduces the
probability of destabilisation of receding liquid rims (Jackiw & Ashgriz,
2022) and production of fine drops (Néel & Villermaux, 2018); but our
following results show that many interesting breakup mechanisms can already be
captured with this choice of $L_{\rm sig}$, which we will present further
below.

(a)

(b)

(c)

Figure 10: Ensemble-averaged instantaneous size distribution functions (a),
and probability distribution functions of axial (b) and radial (c) speed of
aerobreakup fragments at various grid resolution and signature levels. For all
test cases, $We=15$ and $Oh=10^{-3}$.

We first show in fig. 10 the time evolution of the instantaneous distributions
of the size, axial and radial speed distributions of the fragments produced
from bag breakup. Figure 10a indicates that immediately after the onset of bag
breakup ($t/\tau=1.96$) only small fragments close to the minimum grid size
are produced, and well-resolved larger fragments satisfying $d\geq
8\Delta_{\rm sig}$ only come into existence as time elapses, and are always
fewer compared with small fragments near the grid size. The shape of the size
distribution function gradually stabilises, and reaches a steady state by
$t/\tau=2.11$. These findings suggest that the smaller and larger fragments
are produced through different physical mechanisms that arise at different
stages of bag breakup, and eventually these fragmentation mechanisms die out
as the bag approaches full disintegration. The remaining rim will then
disintegrate at still later times, whose investigation we leave for future
work.

Figures 10b and 10c show the instantaneous distribution of fragment axial and
radial speeds $u_{x}$ and $u_{r}\equiv\sqrt{u_{y}^{2}+u_{z}^{2}}$ as functions
of their sizes. It can be seen that the speed of well-resolved fragments
satisfying $d\geq 8\Delta_{\rm sig}$ remains close to a constant value without
significant variations. While the statistics of smaller fragments with $d\leq
8\Delta_{\rm sig}$ are not fully numerically converged, they do show
considerably larger variation around the binned average value in qualitative
agreement with the experimental results of Guildenbecher et al. (2017).
Interestingly, we observe peaks around $d/R_{0}=10^{-2}$ in the distributions
of both $u_{x}$ and $u_{r}$, whose location does not appear to vary with time.
Moreover, fig. 10b suggests that the average axial speed $u_{x}$ of smaller
fragments with $d\leq 8\Delta_{\rm sig}$ increases over time, whereas the
radial speed $u_{r}$ does not show similar increasing trend in fig. 10c. This
is most likely because the smaller fragments are generated earlier and
therefore are exposed to the airflow for much longer periods of time compared
with larger fragments; together with their smaller mass, this means that they
are much more easily accelerated by the axial airflow, hence the continuous
increase in their $u_{x}$ values. On the other hand, $u_{r}$ does not increase
significantly over time, likely because the airflow in the wake region does
not have a large radial velocity component that can accelerate bag fragments
as they migrate downstream.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 11: Snapshots showing the non-local breakup of a long ligament into
multiple fragments during bag film fragmentation with $We=15$ and
$Oh=10^{-3}$. The red boxes shows formation of a single fragment through non-
local end-pinching and its subsequent oscillation, and the blue boxes show the
formation of two fragments through a local breakup event and their subsequent
coalescence.

We will hereafter discuss qualitatively several mechanisms through which the
bag film undergoes fragmentation and form small droplets, which can be
identified by inspecting typical simulation snapshots taken from our
$L=14,\,L_{\rm sig}=13$ simulations. Firstly, fig. 11 shows the non-local
breakup of a stretched long ligament neighbouring two enlarging holes. It can
be seen from fig. 11a that significant cross-sectional diameter variations
have developed on the ligament before the onset of its breakup, which can be
viewed as the result of the nonlinear development of the RP instability (Pal
et al., 2021). Afterwards, the ligament shrinks to form sharp tips and then
breaks up on multiple sites, as shown in fig. 11b, and forms a primary drop
which continues to undergo periodic prolate-oblate shape oscillations
resembling droplets formed out of breaking Rayleigh jets (Hu et al., 2021), as
highlighted in the red boxes in Figs. 11b-11h. This is because the pinch-off
of the stretching ligament induces an inner velocity field within the
detaching droplet that drives it in the oblate direction (see e.g. fig. 7b in
(Hu et al., 2021)), matching the perturbation shape of the second Rayleigh
mode, which then excites oscillation modulated by capillary effects. In the
meantime, the other parts of the ligament do not pinch off to form a series of
fragments at once, but first break up into several elongated structures, and
then form large primary and small satellite drops via the well-known end-
pinching mechanism (Castrejón-Pita et al., 2012; Pal et al., 2021), which
under certain circumstances might coalesce and form a larger fragment as
highlighted in the blue boxes, resembling the ‘immediate satellite merge’
mechanism discussed by Vassallo & Ashgriz (1991).

Figure 12 first shows another type of short ligament breakup and its eventual
evolution into a single droplet, as highlighted in the blue boxes. Compared
with the breakup of long ligaments demonstrated in fig. 11, this type of short
ligament breakup bears stronger resemblance to the breakup phenomena of liquid
bridges studied by Agbaglah (2021), where fragments produced from the same
liquid bridge do not feature significant size variations. This is most likely
because the initial holes are placed very close to each other in Agbaglah
(2021), and the receding rim does not have enough time to grow in size and
momentum before their impact.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 12: Snapshots showing the detachment of a liquid node from ligament
webs (highlighted in the red box) and the evolution of a short ligament into a
single drop (highlighted in the blue box) during bag film fragmentation with
$We=15$ and $Oh=10^{-3}$.

Another type of fragmentation mechanism can also be identified in fig. 12; as
highlighted in the red boxes, three adjacent holes have merged with each
other, and their three bordering rims converge on a common ‘node’ as they are
stretched, which is also observed in the breakup of ligament webs formed on
Savart sheets by Lhuissier & Villermaux (2013). Compared with the ligament
pinch-off mechanism discussed earlier, the surface evolution of this ‘node’
shows much more complicated corrugation patterns as the rims it was connected
to are gradually detached, and therefore is not dominated by the second
Rayleigh mode alone. The ‘node’ drop that eventually forms in this case also
has a much larger size compared with its counterparts formed from ligament
pinch-off events. However, different from an earlier study by Vledouts et al.
(2016), our choice of $L=14$ and $L_{\rm sig}=13$ does not produce holes on
the bag film with as high a number density, and therefore we are not able to
directly measure the wavelength of the RT instability responsible for the film
fragmentation from the average distance between the centres of adjacent holes.
The relatively smaller number density of holes formed also means that the
larger ‘node’ drops formed due to the merging of three or more adjacent and
similarly sized holes are relatively rare compared with generally smaller
fragments formed from ligament breakup, which require the collision of only
two adjacent liquid rims. Furthermore, no less than three holes should fully
expand and arrive at the same region on the bag film where the node is
located, and each of the connecting rims need to break off successively before
the node can be treated as a seperate fragment by the fragment counting
algorithm. These factors help explain the tail of the size distribution
function of aerobreakup fragments taking shape at much later time as we
observed in fig. 10a; namely, that larger droplets are relatively few, and
produced at generally later times during the fragmentation process.

(a)

(b)

(c)

(d)

Figure 13: Snapshots showing the evolution of ‘fingering’ liquid lamellae
during bag film fragmentation with $We=15$ and $Oh=10^{-3}$.

We also note that while we often observe the formation, oscillation and
subsequent corrugation development of liquid ligaments after the impact of
receding rims, as shown on the ligament to the left of the red box in fig.12,
destabilisation of such structure due to the RT instability and its subsequent
evolution into fully developed transverse ‘fingers’ and ‘fine drops’ are only
occasionally observed in our current simulations. According to Néel et al.
(2020), these two regimes are separated by a critical local Weber number for
rim collision $We_{c}\equiv\rho_{l}(2v_{\rm TC})^{2}d_{l}/\sigma=66$, where
$d_{l}$ is the rim diameter. Neglecting the curved geometry of the bag film,
liquid mass conservation further yields $We_{c}=8\sqrt{D_{c}/\pi h}$, where
$D_{c}$ and h are respectively the distance between the centre of two
neighbouring holes and the film thickness. Taking $h=3D/2^{L_{\rm
sig}}=3D/2^{13}$, we find that $We_{c}=66$ corresponds to $D_{c}=2.3R_{0}$,
which we expect might be reached for some pairs of sufficiently separated
holes, as the bag diameter $d_{f}$ before the onset of fragmentation typically
approaches $4R_{0}$, as shown in both our numerical results and the
experimental data of Jackiw & Ashgriz (2021) (see e.g. their fig. 30). Figure
13 highlights two such examples in red boxes, where we observe the transverse
growth of the lamella and the growth of finger-shaped corrugations on its
edges; however, before the fingers fully develop and detach as ‘fine’ drops,
holes are observed to form on the thinning lamella, which then expand and
collide with the fingering lamella edges, turning them into isolated breaking
ligaments. Similar phenomena of lamellae rupture and their edges forming
corrugated ligaments can also be observed in fig. 14 of Vledouts et al.
(2016), although in that case the lamellae appear to remain within the plane
of the film surface, and do not experience transverse growth; and VOF breakup
may play a role in the present examples. Nevertheless, the liquid ligaments
found in the current simulations have already displayed a variety of well-
documented physical phenomena that collectively contribute to the large span
of fragment sizes found in our fragment distribution functions.

(a)

(b)

(c)

(d)

Figure 14: Snapshots showing the receding liquid rim destabilisation during
bag film fragmentation with $We=15$ and $Oh=10^{-3}$. The sites where the rim
is detached from its base is highlighted in red boxes.

Lastly, we observed a few examples showing the destabilisation of receding
liquid rims, as demonstrated in fig. 14. These destabilisation phenomena are
absent in the recent work of Agbaglah (2021) where holes expand over a flat
liquid film, and are therefore most likely linked with the influence of
centrifugal acceleration caused by the curved bag film (Lhuissier &
Villermaux, 2012; Jackiw & Ashgriz, 2022). While we are not yet able to
measure the wavelength or the linear growth rate of the instability, and
therefore have not identified which type of hydrodynamic instability is
involved here; we observe in fig. 14d regular-spaced holes highlighted in the
red boxes forming at the foot of the rim bordering the hole on the left. These
are not attributable to VOF breakup because they are not observed to appear
elsewhere on the bag at the same time; while the rim bordering the larger hole
on the right is seen to develop regular corrugation patterns, which might be
an indication that the receding liquid rim is experiencing the RP instability.
Any further development of the instability is interrupted by the eventual
collision between adjacent rims (not shown in fig. 14). The readers are
referred to Jackiw & Ashgriz (2022) for a more comprehensive discussion on rim
destabilisation. There a few candidates are proposed, including the RT
instability mechanism governing the ‘fingering’ behaviour on bursting surface
bubbles (Lhuissier & Villermaux, 2012); and it is concluded that the
centrifugal acceleration does not govern the instability of the rim directly,
but instead regulates the thickness of the rim via a local-Bond-number
criterion (Wang et al., 2018), with the rim in turn susceptible to the onset
of the RP instability. With the present methodology, these rim instabilities
can be investigated in more detail with higher signature levels $L_{\rm sig}$,
and concomitantly higher $L$. Given the large computational expense of such
simulations, it is not feasible to include such an analysis in the present
study.

###  4.3 Behaviour of bag fragments

In the following, we move away from considering the dynamics and numerical
characteristics of the production mechanisms of fragments to examine instead
those of the fragments themselves. To provide further insight into the
evolution of individual fragments rather than their collective behaviour, we
utilise the droplet tracking algorithm proposed by Chan et al. (2021a) in
post-processing to reconstruct their breakup lineage. This toolbox assumes
breakup and coalescing events to be binary (i.e. at most two parent droplets
may collide and form one large child droplet, or two child droplets may be
produced from the breakup of one parent droplet in a single breakup/coalescing
event), and is capable of identifying all coalesce/breakup events and
differentiate between the new drops produced from these events and those drops
which do not undergo such changes. It requires only the instantaneous fragment
size, location and velocity output from the simulation at given time
intervals, instead of knowledge of the entire flow field at successive
simulation time steps, and therefore incur only limited computational cost
(Chan et al., 2021a).

(a)

(b)

Figure 15: (a): The lifetime of the parent fragments $T_{p}$ as a function of
their diameter $d_{p}$, where $We=15$ and $Oh=10^{-3}$. The bin-averaged
results are shown in grey hollow squares, and the original data are shown as
solid dots, whose colour represents the value of the child/parent diameter
ratio. (b): Velocity difference between parent and child fragments $\Delta u$
as a function of the child/parent diameter ratio $d_{c}/d_{p}$ is shown in the
main plot, whereas the inset plots $\Delta u$ as a function of the diameter of
child fragments $d_{c}$.

As the fragments produced from bag breakup are much smaller compared with the
parent drop, and therefore have a much smaller Weber number, it is unlikely
that they will undergo another bag breakup event. However, they may still
experience secondary breakup to form smaller fragments as they evolve over
time. It is therefore of interest to determine the lifetime of breaking parent
fragments $T_{p}$ using the toolbox of Chan et al. (2021a), defined as the
interval between their birth and death in two successive breakup events
(Rivière et al., 2022). Figure 15a shows $T_{p}$ as a function of the diameter
$d_{p}$ of parent fragments, with the bin-averaged values of $T_{p}$ shown in
grey squares. The solid dots plotted in the background represent recorded
individual breakup events, and are colour-coded by the logarithm value of the
minimum child/parent diameter ratio $d_{c}/d_{p}$, highlighting a broad
distribution of parent fragment lifetime. For comparison, the characteristic
capillary time of fragments $t_{\rm cap}$ is also plotted in fig. 15a as a
function of $d_{p}$, which is defined as follows:

$t_{\rm cap}(d_{p})=\sqrt{\frac{\rho_{l}d_{p}}{2\sigma}}.$ (9)

It is noted that breakup events satisfying $d_{p}/R_{0}\approx 2$ in fig. 15a
correspond to the primary breakup of the bag film still attached to the
remnant rim, while those satisfying $d_{p}/R_{0}\leq 0.6$ represent the
secondary breakup of bag fragments. It can first be seen from the scattered
original data that most of the bag film fragments that undergo a secondary
breakup fall within the range of $d_{p}\geq 8\Delta_{\rm sig}$, which
correspond to the ‘well resolved’ fragments discussed in §4.1. Furthermore,
the lifetime of fragments $T_{p}$ satisfying $d_{p}\leq 0.1R_{0}$ shows a
dependency on $r_{p}$ that roughly scales with the characteristic capillary
time, but this trend breaks down for even larger fragments with $d_{p}\geq
0.1R_{0}$. It is noted that the capillary time defined in Eq. (9) is
proportional to the oscillation period of droplet spherical harmonic modes
(Eq. (12)). Therefore, the scaling of $T_{p}$ with $t_{\rm cap}$ for
$d_{p}\leq 0.1R_{0}$ may suggest that the fragmentation of these fragments is
primarily due to large-amplitude nonlinear oscillations which can trigger a
capillary breakup (Lalanne et al., 2019). This also explains the large
scatters we observe in the lifetime of parent fragments within this size
range, as when nonlinearity becomes dominant, the surface oscillations cannot
be represented by a single mode, and different modes of perturbation with
different oscillation periods might trigger breakup depending on specific
fragments. As for even larger fragments with $d_{p}\geq 0.1R_{0}$, their
lifetime $T_{p}$ appears to scatter around an average value of $10t_{\rm
cap}(h_{c})$, where $t_{\rm cap}(h_{c})$ is the characteristic capillary time
based on the critical film thickness $h_{c}=3D/2^{L_{\rm sig}}$. Here
$10t_{\rm cap}(h_{c})$ is an estimation of the inertial timescale leading to
the capillary breakup of stretching liquid ligaments, which are formed due to
hole collision (Agbaglah, 2021). This suggests that this type of non-local
breakup is only dependent on the topological evolution of the stretching
liquid ligament, rather than that of the entire parent drop from which the
child fragments are torn off; but this remains to be verified in future work.

We were also able to compute the magnitude of the velocity differences $\Delta
u$ between fragment parents and their children at two successive instants when
the fragment statistics are collected, and plot them in fig. 15b as a function
of the ratio between the child and parent diameter $d_{c}/d_{p}$. It is found
that for breakup events where a small child/parent size ratio
$(d_{c}/d_{p}\leq 0.22)$ are detected, the velocity difference $\Delta u$
appear to show little dependence on $d_{c}/d_{p}$, despite significant
scatter. Based on our findings in fig. 15a, these breakup events with
$d_{c}/d_{p}\leq 0.22$ mostly feature small children with large parents. On
the other hand, breakup events satisfying $d_{c}/d_{p}\geq 0.22$ are dominated
by small parents and children, and their $\Delta u$ decreases with increasing
values of $d_{c}/d_{p}$, roughly following a power-law scaling $\Delta
u\propto(d_{c}/d_{p})^{-1.2}$. We further note that the $\Delta u$ values of
breakup events satisfying $d_{c}/d_{p}\leq 0.22$ roughly coincides with the
inviscid Taylor-Culick velocity $v_{\rm TC}$, with an estimation for the bag
film thickness $h$ based on the signature level $L_{\rm sig}$:

$v_{\rm TC}\equiv\sqrt{\frac{2\sigma}{\rho_{l}h}}=\sqrt{\frac{2^{L_{\rm
sig}+1}\sigma}{3\rho_{l}D}}\approx 0.17U_{0}.$ (10)

This agrees with the recent confirmation by Néel & Deike (2022) that the speed
of film drops produced from bubble-bursting can be estimated by $v_{\rm TC}$
to an order of magnitude. An explanation for this approximate agreement is as
follows. Prior and up to collision between adjacent hole rims, each rim
travels at the Taylor-Culick velocity. The colliding rims of these holes then
form liquid ligaments, which exhibit an axial stretching rate comparable to
the pre-collision rim speed (i.e. the Taylor-Culick velocity). This stretching
rate in turn sets the relative speed of sufficiently small fragments ejected
from the parent ligament. The parent ligaments may constitute part of the
parent drop, or may themselves have separated from it. Therefore, for a given
child droplet size produced by this mechanism, the ratio $d_{c}/d_{p}$ may see
considerable variation. Consequently, we observe the large range of
$d_{c}/d_{p}\leq 10^{-2}$ where the parent/child velocity difference remains
close to $v_{\rm TC}$. The inset of fig. 15b further shows the velocity
difference in breakup events as a function of the child diameter $d_{c}$
alone; and this time we find that it is the small fragments satisfying
$d_{c}\leq 8\Delta_{\rm sig}$ that appear to approach $v_{\rm TC}$, agreeing
with our analysis that it is the fragments produced from liquid ligament
breakup that are represented within the range of $d_{c}/d_{p}\leq 0.22$. Note
however that in the inset of fig. 15b, $\Delta u$ continues to increase beyond
$v_{\rm TC}$ with decreasing $d_{c}$, and does not scatter around $v_{\rm TC}$
as in the main plot. This is most likely because the droplet tracking
algorithm has a finite fragment detection frequency (Chan et al., 2021a), and
the smallest fragments may be generated and accelerated by the ambient airflow
between two successive instants when this algorithm is called for their
detection. The nominal velocity difference $\Delta u$ for the smallest
fragments therefore includes contributions from both the initial ejection
velocity (primarily in the radial direction) and the increment due to airflow
acceleration (primarily in the axial direction), which causes $\Delta u$ to
increase steadily beyond $u_{\rm TC}$. As these smallest fragments are
scattered in different averaging bins according to the child/parent diameter
ratio of the breakup events, the contribution of airflow acceleration to
$\Delta u$ becomes less obvious.

Based on our findings in fig. 15b, we introduce the following criteria for
determining ‘non-local’ breakup and coalescing events, respectively:

$\max\left(\frac{d_{c,i}}{d_{p}}\right)\leq
0.22,\quad\max\left(\frac{d_{c}}{d_{p,i}}\right)\leq 4.64,\quad i=1,2.$ (11)

Otherwise we term the breakup or coalescing events ‘local’. Here the critical
child/parent diameter ratio of 0.22 for fragment breakup (and analogously 4.64
for coalescence) separates the two breakup regimes found in fig. 15b, where
the velocity difference $\Delta u$ either scatters around $v_{\rm TC}$ or
scales with $d_{c}/d_{p}$. Note that the term ‘local’ here does not mean the
parent and child fragments are close to each other in terms of their locations
in the physical space, which all such events satisfy; but rather in the sense
of the parent and its two children being close in their respective sizes (Chan
et al., 2021b).

(a)

(b)

Figure 16: Ensemble-averaged evolution of number of breakup (a) and coalesce
(b) events during the breakup of bag films produced from an initial droplet
with $We=15$ and $Oh=10^{-3}$. All ensemble realisations are run at
$L=14,\,L_{\rm sig}=13$.

We further plot in fig. 16 the ensemble-averaged number density of breakup and
coalescing events detected during the bag fragmentation period. It is found
that both breakup and coalescing events occur most frequently around
$t/\tau=2.10$. This is most likely when the receding rims fully absorb the bag
film and collide with each other, which then triggers a series of corrugated
ligament breakup and fragment coalescing events. After this the breakup and
coalescing behaviour become less frequent, as the corrugated ligaments
gradually disintegrate without liquid mass input from the bag film. While
there exists ‘multistep’ breakup in our aerobreakup simulation outputs, our
results in fig. 16 suggest that the fragmentation process involved in the
aerobreakup problem cannot be well described by a breakup cascade model
(Garrett et al., 2000; Chan et al., 2021b), as non-local breakup and
coalescing events producing children with sizes drastically different from
their parents are found to dominate this problem, which is different from the
entrained air bubble breakup scenario in breaking wave studies (Garrett et
al., 2000; Deane & Stokes, 2002; Deike et al., 2016; Chan et al., 2021b;
Mostert et al., 2022) where the prevalence of local breakup events leads to a
well-defined bubble-mass flux supporting breakup cascade models. It is also
noted that breakup events occur much more frequently than coalescing events
for bag films, which is expected as the latter requires two adjacent fragments
to cross paths at the same time, which only happen for a small portion of
neighbouring fragments with specific initial position and velocity
configurations.

(a)

(b)

Figure 17: Oscillatory behaviour of very small aerobreakup fragments produced
from an initial droplet with $We=15$ and $Oh=10^{-3}$. Left: surface energy
evolution of individual fragments (blue curves) with their steady-state
surface energy plotted (dashed lines) for reference; right: frequency of the
dominant fragment oscillation mode as a function of the fragment radius.

Apart from identifying all breakup and coalescing events in the spray formed
due to aerobreakup, the toolbox developed by Chan et al. (2021a) also enables
us to track the evolution of properties of individual fragments during their
lifetime. For example, fig. 17a shows the evolution of surface energy $E_{s}$
of a portion of individual small fragments. Decaying oscillation patterns are
clearly visible in fig. 17a, with their oscillation frequency generally
increasing with decreasing fragment radius $r$ (hence decreasing steady-state
surface energy values). We further extract the frequency of the dominant
oscillation mode of small fragments at two different grid resolution levels
($L=14,\,L_{\rm sig}=13$ and $L=13,\,L_{\rm sig}=11$), and plot them against
the fragment radius in fig. 17b, where an excellent agreement is found for
small fragments within the diameter range of $0.02R_{0}\leq d\leq 0.2R_{0}$
with the theoretical predictions of Prosperetti (1980) for the second Rayleigh
mode, which is given as follows for an inviscid droplet with density
$\rho_{l}$ and radius $R^{*}$ at equilibrium:

$\omega_{n,0}=\sqrt{(n-1)n(n+2)\frac{\sigma}{\rho_{l}R^{*}}},$ (12)

where $n$ is the spherical harmonic mode number at which the interface of the
droplet is perturbed. The second Rayleigh mode corresponds with $n=2$, and
this mode number is associated with the oblate-prolate shape perturbations
which we observed in fig. 11. Both the viscous and inviscid theoretical model
of Prosperetti (1980) are plotted in fig. 17, which almost completely overlap
except for small fragments below $8\Delta_{13}$, where the viscous model shows
a slightly better match with the numerical results. This is because the $Oh$
value of $10^{-3}$ at which simulations discussed in this section are run is
very low, and viscous effects become non-trivial only for very small
fragments. Furthermore, for results at both resolution levels shown in fig.
17b, the agreement between numerical and theoretical results reaches into
their corresponding range of small fragments with $d\leq 8\Delta_{\rm sig}$,
although the fragment size and velocity distributions within this range have
not reached grid convergence. Overall, these results demonstrate that
fragments satisfying $d\leq 8\Delta_{\rm sig}$ are found to be governed by
well-documented physical mechanisms, e.g. rim retraction at the Taylor-Culick
velocity and the Rayleigh oscillation theory (Prosperetti, 1980), even though
full grid independence in terms of fragment size and velocity statistics is
still to be established. This suggests that it is the droplet production
mechanism through ligament fragmentation which remains somewhat grid-
dependent, while the droplets resulting from these production events are
themselves well-resolved. Nevertheless, we highlight that the film
fragmentation process is well-resolved in our numerical simulations with the
aid of the MD algorithm (Chirco et al., 2022).

###  4.4 Viscous effects on bag breakup

Having provided an overview of the physical mechanisms governing bag breakup
and the subsequent evolution patterns of the fragments at a specific $Oh$
value of $10^{-3}$ in §4.2, in this section we further vary the $Oh$ value and
examine its influence on the bag breakup phenomena. This has not been examined
in depth in currently available aerobreakup studies as most research efforts
have been carried out in the limit of very low $Oh$ values (Guildenbecher et
al., 2017; Jackiw & Ashgriz, 2022; Kant et al., 2022).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 18: Simulation snapshots showing the bag breakup process at different
$Oh$ values ($10^{-4}$, $10^{-3}$, $10^{-2}$ and $5\times 10^{-2}$ from the
top to the bottom row), where $We$ is fixed as 15. For all cases, $L=14$ and
$L_{\rm sig}=13$.

We first provide in fig. 18 simulation snapshots showing the bag breakup
process for $Oh=10^{-4}$, $10^{-3}$, $10^{-2}$ and $5\times 10^{-2}$ for a
qualitative analysis, with $We=15$ for all cases. It first be seen that as
$Oh$ increases, the bag becomes more ‘flattened’ and its surface area becomes
smaller, and correspondingly the surrounding rim around the bag becomes more
prominent. This implies that the inviscid model proposed by Jackiw & Ashgriz
(2021) predicting the volume of the bag film and rim may need to be extended
for a generalisation to the highly viscous regime. Furthermore, it is observed
that the ligament breakup behaviour changes significantly as the $Oh$ value
increases. While at $Oh=10^{-4}$ the receding liquid rims generate capillary
waves propagating through the entire bag (Savva & Bush, 2009), and undergo
destabilisation patterns similar to what we observed in fig. 14, these are not
found at higher $Oh$ values, which suggests that these phenomena are highly
sensitive to viscous damping effects, and their contribution to fragment
statistics becomes negligible with increasing $Oh$. At $Oh=10^{-3}$, liquid
ligaments typically show long periods of radial oscillation after being formed
out of colliding hole rims, and then break up to form ‘primary’ and
‘satellite’ drops that differ significantly in their sizes. When $Oh$
increases to $10^{-2}$, it is found that far fewer satellite drops are
produced from ligament breakup, and the ‘end-pinching’ breakup mechanism comes
to dominate as the ligaments now tend to break up on one end repeatedly and
form small drops. This may be because the ligaments can be stretched longer
and thinner with a higher $Oh$, and the smaller ligament radius impedes the
formation of satellite drops. According to Vassallo & Ashgriz (1991), smaller
radius of ligaments induces a larger pressure difference that pushes the free
end back in the axial direction much quicker, hence preventing capillary
pinch-off in the radial direction that produces the satellite drops. Further
increasing $Oh$ to $5\times 10^{-2}$ causes the ligaments to be stretched even
thinner and produce smaller fragments once they break up, which is because
increased viscosity smooths out the variation of the axial velocity along the
ligament that drive the pinch-off events (Hu et al., 2021). Another side
effect appearing at $Oh=5\times 10^{-2}$ is that fewer node fragments are
observed to form, which is because the decreased bag area leaves smaller room
for generation and mutual collision of more than three holes which produce
node fragments. Savva & Bush (2009) suggest that at even higher $Oh$ values
the liquid rims will disappear, and the thickness of the entire bag film will
correspondingly increase as the holes enlarge, although we do not reach this
limit in our current numerical simulations.

(a)

(b)

(c)

Figure 19: (a): Time evolution of the total number of film fragments after the
onset of bag breakup for one ensemble realisation with different $Oh$ values.
(b)(c): Time- and ensemble-averaged size (b) and speed (c) probability
distribution functions of aerobreakup fragments with $We=15$ and various $Oh$
values.

Figure 19 further shows the evolution of the number of fragments satisfying
$d\geq\Delta_{\rm sig}$, and the time- and ensemble-averaged size and speed
distribution functions for the three specific $Oh$ values selected in fig. 18.
Figure 19a shows that as the $Oh$ value increases, the number of fragments
$N_{\rm frag}$ reached when the bag fully disintegrates decreases, which is
because the total area of the bag film decreases, leaving less amount of
liquid that feed the film breakup process. While $N_{\rm frag}$ generally
increases over time despite small-scale local oscillations, which most likely
arise from relatively rare coalescing events, long periods of time where
$N_{\rm frag}$ remains nearly constant can be clearly seen for $Oh=5\times
10^{-2}$, as the highly-viscous liquid ligaments can now be sustained for much
longer under stretching, and it is their intermittent breakup events that
contribute to the isolated sharp growth events in $N_{\rm frag}$. Figure 19b
shows that the fragment size distribution functions for $10^{-4}\leq Oh\leq
10^{-2}$ remain very close to each other. When $Oh$ is further increased to
$5\times 10^{-2}$, it is observed that the size distribution function for
$d\geq\Delta_{\rm sig}$ becomes much more convex-shaped, with more
‘intermediate’ fragments produced within the range of $\Delta_{\rm sig}\leq
d\leq 8\Delta_{\rm sig}$, and fewer large fragments with $d\geq 8\Delta_{\rm
sig}$. This arises from the coupled effects of reduction of bag film area
(hence smaller chance for formation of ‘node’ fragments, as observed in fig.
18l) and breakup of long viscous ligaments into smaller fragments, and implies
a reduction in the average fragment size which will be discussed in more
details further below. Finally, fig. 19c suggests that increasing the $Oh$
value causes the fragment speed to become more evenly-distributed across
different fragment sizes, which is probably due to the combined effects of
viscous damping of the internal axial velocity distributions of pre-breakup
ligaments (Hu et al., 2021), and the airflow becoming much less turbulent as
its viscosity also increases under a fixed viscosity ratio $\mu^{*}$.

(a)

(b)

Figure 20: The breakup onset time $t_{b}$ (a) and the instantaneous average
diameter $\bar{d}$ of fragments satisfying $d\geq 8\Delta_{\rm sig}$ (b) as
functions of the $Oh$ values, with an exponential model fit for (a). The non-
dimensionalised experimental data of Kant et al. (2022) are included in (b)
for comparison. Squares mean the results have been ensemble-averaged over
three individual realisations, and crosses mean data from only one realisation
is available. For all cases $We=15$.

We now analyse the influence of $Oh$ on the ensemble-averaged breakup time
$\bar{t}_{b}$ (defined as the time when the first hole is generated on the
film by the MD algorithm) and the ensemble-averaged instantaneous diameter of
fragments $\bar{d}_{\rm frag}$ in fig. 20. It should be noted that due to the
significant runtime required on supercomputers, the results at a few $Oh$
values in fig. 20 are obtained from only one realisation instead of being
ensemble-averaged, and these results are differentiated from the others by a
cross mark. Nevertheless, from fig. 20a, it is seen that $\bar{t}_{b}$ first
remains almost independent of $Oh$ as the latter increases from $10^{-4}$ to
$10^{-2}$, which is likely because the wake region behind the droplet is fully
turbulent, and the thinning process of the bag before the onset of its breakup
is mainly influenced by capillary and inertial effects, as discussed in §3.2.
When $Oh$ exceeds $10^{-2}$, $\bar{t}_{b}$ is found to increase dramatically
with $Oh$ in an exponential manner as shown by the fitted model, which may be
the consequence of both the transition of the wake region from a fully
turbulent to a laminar status (hence smaller fore-aft pressure difference on
the deforming drop that pushes out the bag), and the bag thinning process
coming under the domination of a capillary-viscous balance. The hypothesis of
the influence of wake region on the growth of $t_{b}$ is further supported by
examining the freestream Reynolds number $Re$:

$Re\equiv\frac{\rho_{a}U_{0}d_{0}}{\mu_{a}}=\frac{\mu^{*}\sqrt{We}}{\sqrt{\rho^{*}}Oh}=\frac{7.381}{Oh},$
(13)

where the critical $Oh$ value of $10^{-2}$ corresponds to $Re=7.38\times
10^{2}$, which agrees with the order of magnitude of previously-reported $Re$
values at which the wake regions behind a sphere transits to turbulence and
vortex shedding is initiated (Rodriguez et al., 2011). Overall, this sharp
increase of the breakup time $\bar{t}_{b}$ with increasing $Oh$ beyond
$10^{-2}$ agrees with the early findings that the $Oh$ values do not have
significant influence over breakup regimes when they are below 0.1 (Hsiang &
Faeth, 1992).

Finally, we compute the instantaneous average diameter $\bar{d}_{\rm frag}$ at
the end of bag breakup. The averaging is completed for fragments satisfying
$d\geq 8\Delta_{\rm sig}$ over different individual realisations with the same
$Oh$ value, which enables us to acquire sufficient amounts of numerically
converged statistics to produce meaningful results. In fig. 20b, $\bar{d}_{\rm
frag}$ shows a non-monotonic dependence on the film Ohnesorge number, defined
as $Oh_{\rm film}\equiv\mu_{l}/\sqrt{\rho_{l}\sigma h_{f}}$. While the
thinning process of bag films continues as they undergo fragmentation, we
select a constant characteristic film thickness value $h_{f}=D/2^{L_{\rm
sig}}$ so that our computed average fragment diameters match the order of
magnitude of the fragment statistics produced from the breakup of thin films
studied by Kant et al. (2022), which are also included in fig. 20b for
comparison. Namely, as $Oh_{\rm film}$ increases from $1.65\times 10^{-3}$ to
0.826, $\bar{d}_{\rm frag}$ first remains close to $25h_{f}$, followed by an
abrupt decrease as it approaches $Oh_{\rm film}=0.4$. This non-monotonic
dependency on $Oh$ is also observed in the results of Kant et al. (2022),
where an initial increase of $\bar{d}_{\rm frag}$ to $29.2h_{f}$ is followed
by an abrupt decrease to $14.2h_{f}$ when $Oh_{\rm film}$ increases beyond
unity. Based on our analysis of fig. 18 and fig. 19, we ascribe the abrupt
decrease of $\bar{d}_{\rm frag}$ with increasing $Oh$ values to the formation
of much fewer large ‘node’ fragments due to the decrease of bag area, and the
breakup of ligaments that are stretched much thinner under high viscosity.

##  5 Conclusions

We have presented in this study the results of both axisymmetric and three-
dimensional numerical simulations of droplet aerobreakup. For the axisymmetric
simulations, our results were validated by a good agreement with the constant
spanwise growth rate of drops predicted by Jackiw & Ashgriz (2021), and we
were able to explain deviation from this model by the interaction between the
drop surface and the wake vortices. We were also able to look into the
thinning of bag films before the onset of bag breakup, and found that at small
$Oh$ values capillary effects will cause the thinning rate to exceed that
predicted by Villermaux & Bossa (2009).

For the three-dimensional aerobreakup simulations, we utilised the MD
algorithm (Chirco et al., 2022) for artificial perforation of thin films,
which enabled us to minimise pollution of fragment statistics by spurious
numerical breakup, and establish grid convergence of fragment statistics for
aerobreakup studies for the first time. Afterwards, we analysed the output
fragment statistics, and were able to reconstruct the breakup lineage and
evolution of individual fragment properties using the postprocessing toolbox
proposed in Chan et al. (2021a). It is found that smaller fragments with their
diameters satisfying $d\leq 8\Delta_{\rm sig}$ are produced from non-local
liquid ligament breakup and tend to undergo decaying surface oscillations
dominated by the second Rayleigh mode, with their ejection velocity set by the
colliding liquid rims receding at the Taylor-Culick velocity; while larger
fragments satisfying $8\Delta_{\rm sig}\leq d\leq 10^{-1}R_{0}$ are formed
from breakup of ‘nodes’ bordering two or more holes, and undergo secondary
local breakup events due to large-amplitude nonlinear oscillations.
Destabilisation of receding rims is also found in some individual
realisations, although they do not contribute significantly to fragment
production under current simulation configurations.

We find in particular that numerical convergence of bag-breakup problems
feature subtle numerical convergence properties:

1\. Without the MD algorithm, numerical convergence cannot be achieved for
thin film fragmentation owing to the VOF-breakup phenomenon.

2\. With the MD algorithm, the production of fragments through thin film
fragmentation shows grid convergence for droplet children with diameter
$d>8\Delta_{\rm sig}$.

3\. With or without the MD algorithm, the production of small droplets close
to the resolution limit resulting from ligament fragmentation occurs
independent of VOF- or MD-induced breakup, and numerical convergence for these
production mechanisms is yet to be established in the present study.

4\. However, once they are produced, the subsequent evolution of small
fragments is well-resolved in the present simulations, even for small
fragments approaching the grid resolution.

Finally, we investigated the influence of $Oh$ on bag film breakup, and it is
found that increasing $Oh$ causes the bag area to decrease and the liquid
ligaments to be stretched much thinner, which generally lead to production of
fewer fragments with smaller average diameters.

Overall, these results show the utility of the MD algorithm in improving the
grid convergence behaviour of fragment statistics and helping to recover
previously unresolved fluid physics in two-phase numerical simulations
involving breakup of thin films (Kant et al., 2022), and also shed light on
the effect of viscosity on bag breakup which has not been discussed in detail
in previous aerobreakup studies (Jackiw & Ashgriz, 2022). They also pave the
way for future studies investigating the later development of the remnant rim
and the effects of airphase turbulence and initial perturbations on the
deformation and breakup of droplets, while also serving as a stepping stone
towards a full-scale numerical study of spume drop generation on the air-sea
interface under high wind conditions.

##  6 Declaration of Interests

The authors report no conflict of interest.

##  7 Acknowledgments

The authors would like to thank EPSRC for the computational time made
available on the UK supercomputing facility ARCHER2 via the UK Turbulence
Consortium (EP/R029326/1). Use of the University of Oxford Advanced Research
Computing (ARC) facility is also acknowledged. K. Tang is supported by a
Research Studentship at the University of Oxford.

## References

  * Agbaglah (2021) Agbaglah, G. G. 2021 Breakup of thin liquid sheets through hole–hole and hole–rim merging. Journal of Fluid Mechanics 911.
  * Bourouiba (2021) Bourouiba, L. 2021 The fluid dynamics of disease transmission. Annual Review of Fluid Mechanics 53, 473–508.
  * Brackbill et al. (1992) Brackbill, J. U., Kothe, D. B. & Zemach, C. 1992 A continuum method for modeling surface tension. Journal of Computational Physics 100 (2), 335–354.
  * Castrejón-Pita et al. (2012) Castrejón-Pita, A. A., Castrejón-Pita, J. R. & Hutchings, I. M. 2012 Breakup of liquid filaments. Physical Review Letters 108 (7), 074506.
  * Chan et al. (2021a) Chan, W. H. R., Dodd, M. S., Johnson, P. L. & Moin, P. 2021a Identifying and tracking bubbles and drops in simulations: A toolbox for obtaining sizes, lineages, and breakup and coalescence statistics. Journal of Computational Physics 432, 110156.
  * Chan et al. (2021b) Chan, W. H. R., Johnson, P. L. & Moin, P. 2021b The turbulent bubble break-up cascade. part 1. theoretical developments. Journal of Fluid Mechanics 912.
  * Chirco et al. (2022) Chirco, L., Maarek, J., Popinet, S. & Zaleski, S. 2022 Manifold death: a volume of fluid implementation of controlled topological changes in thin sheets by the signature method. Journal of Computational Physics 467, 111468.
  * Chou & Faeth (1998) Chou, W.-H. & Faeth, G. M. 1998 Temporal properties of secondary drop breakup in the bag breakup regime. International Journal of Multiphase Flow 24, 889–912.
  * Coyajee & Boersma (2009) Coyajee, E. & Boersma, B. J. 2009 Numerical simulation of drop impact on a liquid–liquid interface with a multiple marker front-capturing method. Journal of Computational Physics 228 (12), 4444–4467.
  * Deane & Stokes (2002) Deane, G. B. & Stokes, M. D. 2002 Scale dependence of bubble creation mechanisms in breaking waves. Nature 418 (6900), 839–844.
  * Deike et al. (2016) Deike, L., Melville, W. K. & Popinet, S. 2016 Air entrainment and bubble statistics in breaking waves. Journal of Fluid Mechanics 801, 91–129.
  * Erinin et al. (2019) Erinin, M. A., Wang, S. D., Liu, R., Towle, D., Liu, X. & Duncan, J. H. 2019 Spray generation by a plunging breaker. Geophysical Research Letters 46, 8244–8251.
  * Garrett et al. (2000) Garrett, C., Li, M. & Farmer, D. 2000 The connection between bubble size spectra and energy dissipation rates in the upper ocean. Journal of Physical Oceanography 30 (9), 2163–2171.
  * Gorokhovski & Herrmann (2008) Gorokhovski, M. & Herrmann, M. 2008 Modeling primary atomization. Annual Review of Fluid Mechanics 40, 343–366.
  * Guémas et al. (2015) Guémas, M., Sellier, A. & Pigeonneau, F. 2015 Slow viscous gravity-driven interaction between a bubble and a free surface with unequal surface tensions. Physics of Fluids 27 (4), 043102.
  * Guildenbecher et al. (2017) Guildenbecher, D. R., Gao, J. & Chen, J.and Sojka, P. E. 2017 Characterization of drop aerodynamic fragmentation in the bag and sheet-thinning regimes by crossed-beam, two-view, digital in-line holography. International Journal of Multiphase Flow 94, 107–122.
  * Guildenbecher et al. (2009) Guildenbecher, D. R., López-Rivera, C. & Sojka, P. E. 2009 Secondary atomization. Experiments in Fluids 46, 371–402.
  * Hsiang & Faeth (1992) Hsiang, L.-P. & Faeth, G. M. 1992 Near-limit drop deformation and secondary breakup. Int. J. Muhiphase Flow 18, 635–652.
  * Hsiang & Faeth (1995) Hsiang, L.-P. & Faeth, G. M. 1995 Drop deformation and breakup due to shock wave and steady disturbances. International Journal of Multiphase Flow 21, 545–560.
  * Hu et al. (2021) Hu, L., She, L., Fang, Y., Su, R. & Fu, X. 2021 Deformation characteristics of droplet generated by rayleigh jet breakup. AIP Advances 11 (4), 045310.
  * Jackiw & Ashgriz (2021) Jackiw, I. M. & Ashgriz, N. 2021 On aerodynamic droplet breakup. Journal of Fluid Mechanics 913.
  * Jackiw & Ashgriz (2022) Jackiw, I. M. & Ashgriz, N. 2022 Prediction of the droplet size distribution in aerodynamic droplet breakup. Journal of Fluid Mechanics 940.
  * Jain et al. (2015) Jain, M., Prakash, R. S., Tomar, G. & Ravikrishna, R. V. 2015 Secondary breakup of a drop at moderate weber numbers. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 471 (2177), 20140930.
  * Jain et al. (2019) Jain, S. S., Tyagi, N., Prakash, R. S., Ravikrishna, R. V. & Tomar, G. 2019 Secondary breakup of drops at moderate weber numbers: Effect of density ratio and reynolds number. International Journal of Multiphase Flow 117, 25–41.
  * Jalaal & Mehravaran (2012) Jalaal, M. & Mehravaran, K. 2012 Fragmentation of falling liquid droplets in bag breakup mode. International Journal of Multiphase Flow 47, 115–132.
  * Jalaal & Mehravaran (2014) Jalaal, M. & Mehravaran, K. 2014 Transient growth of droplet instabilities in a stream. Physics of Fluids 26 (1), 012101\.
  * Kant et al. (2022) Kant, P., Pairetti, C., Saade, Y., Popinet, S., Zaleski, S. & Lohse, D. 2022 Bags mediated film atomization in a cough machine. arXiv preprint arXiv:2202.13949 .
  * Kočárková et al. (2013) Kočárková, H., Rouyer, F. & Pigeonneau, F. 2013 Film drainage of viscous liquid on top of bare bubble: Influence of the bond number. Physics of Fluids 25 (2), 022105\.
  * Kundu et al. (2012) Kundu, P.K., Cohen, I.M. & Dowling, D.R. 2012 Fluid Mechanics. Elsevier Science.
  * Kékesi et al. (2014) Kékesi, T., Amberg, G. & Wittberg, L. P. 2014 Drop deformation and breakup. International Journal of Multiphase Flow 66, 1–10.
  * Lalanne et al. (2019) Lalanne, B., Masbernat, O. & Risso, F. 2019 A model for drop and bubble breakup frequency based on turbulence spectra. AIChE Journal 65 (1), 347–359.
  * Lhuissier & Villermaux (2012) Lhuissier, H. & Villermaux, E. 2012 Bursting bubble aerosols. Journal of Fluid Mechanics 696, 5–44.
  * Lhuissier & Villermaux (2013) Lhuissier, H. & Villermaux, E. 2013 ‘effervescent’atomization in two dimensions. Journal of Fluid Mechanics 714, 361–392.
  * Li et al. (2019) Li, Y., Zhang, P. & Kang, N. 2019 Theoretical analysis of rayleigh–taylor instability on a spherical droplet in a gas stream. Applied Mathematical Modelling 67, 634–644.
  * Ling et al. (2017) Ling, Y., Fuster, D., Zaleski, S. & Tryggvason, G. 2017 Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: a numerical closeup. Physical Review Fluids 2 (1), 014005.
  * Ling et al. (2015) Ling, Y., Zaleski, S. & Scardovelli, R. 2015 Multiscale simulation of atomization with small droplets represented by a lagrangian point-particle model. International Journal of Multiphase Flow 76, 122–143.
  * Magnaudet & Mercier (2020) Magnaudet, J. & Mercier, M. J. 2020 Particles, drops, and bubbles moving across sharp interfaces and stratified layers. Annual Review of Fluid Mechanics 52, 61–91.
  * Marcotte & Zaleski (2019) Marcotte, F. & Zaleski, S. 2019 Density contrast matters for drop fragmentation thresholds at low ohnesorge number. Physical Review Fluids 4 (10), 103604.
  * Mehta et al. (2017) Mehta, P., Haj-Ahmad, R., Rasekh, M., Arshad, M. S., Smith, A., van der Merwe, S. M., Li, X., Chang, M.-W. & Ahmad, Z. 2017 Pharmaceutical and biomaterial engineering via electrohydrodynamic atomization technologies. Drug Discovery Today 22 (1), 157–165.
  * Mostert & Deike (2020) Mostert, W. & Deike, L. 2020 Inertial energy dissipation in shallow-water breaking waves. Journal of Fluid Mechanics 890.
  * Mostert et al. (2022) Mostert, W., Popinet, S. & Deike, L. 2022 High-resolution direct simulation of deep water breaking waves: transition to turbulence, bubbles and droplets production. Journal of Fluid Mechanics 942.
  * Néel & Deike (2022) Néel, B. & Deike, L. 2022 Velocity and size quantification of drops in single and collective bursting bubbles experiments. Physical Review Fluids 7 (10), 103603.
  * Néel et al. (2020) Néel, B., Lhuissier, H. & Villermaux, E. 2020 ‘fines’ from the collision of liquid rims. Journal of Fluid Mechanics 893.
  * Néel & Villermaux (2018) Néel, B. & Villermaux, E. 2018 The spontaneous puncture of thick liquid films. Journal of Fluid Mechanics 838, 192–221.
  * Nicholls & Ranger (1969) Nicholls, J. A. & Ranger, A. A. 1969 Aerodynamic shattering of liquid drops. AIAA Journal 7 (2), 285–290.
  * Obenauf & Sojka (2021) Obenauf, D. G. & Sojka, P. E. 2021 Theoretical deformation modeling and drop size prediction in the multimode breakup regime. Physics of Fluids 33 (9), 092113.
  * Opfer et al. (2014) Opfer, L., Roisman, I. V., Venzmer, J., Klostermann, M. & Tropea, C. 2014 Droplet-air collision dynamics: Evolution of the film thickness. Physical Review E 89 (1), 013023\.
  * Pairetti et al. (2018) Pairetti, C., Popinet, S., Damián, S., Nigro, N. & Zaleski, S. 2018 Bag mode breakup simulations of a single liquid droplet. In 6th European Conference on Computational Mechanics.
  * Pal et al. (2021) Pal, S., Crialesi-Esposito, M., Fuster, D. & Zaleski, S. 2021 Statistics of drops generated from ensembles of randomly corrugated ligaments. arXiv preprint arXiv:2106.16192 .
  * Pigeonneau & Sellier (2011) Pigeonneau, F. & Sellier, A. 2011 Low-reynolds-number gravity-driven migration and deformation of bubbles near a free surface. Physics of Fluids 23 (9), 092102.
  * Popinet (2009) Popinet, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. Journal of Computational Physics 228 (16), 5838–5866.
  * Popinet (2018) Popinet, S. 2018 Numerical models of surface tension. Annual Review of Fluid Mechanics 50, 49–75.
  * Popinet (2019) Popinet, S. 2019 Basilisk flow solver and pde library. Available at: http://basilisk.fr.
  * Prosperetti (1980) Prosperetti, A. 1980 Free oscillations of drops and bubbles: the initial-value problem. Journal of Fluid Mechanics 100 (2), 333–347.
  * Radhakrishna et al. (2021) Radhakrishna, V., Shang, W., Yao, L., Chen, J. & Sojka, P. E. 2021 Experimental characterization of secondary atomization at high ohnesorge numbers. International Journal of Multiphase Flow 138, 103591.
  * Rivière et al. (2022) Rivière, A., Ruth, D. J., Mostert, W., Deike, L. & Perrard, S. 2022 Capillary driven fragmentation of large gas bubbles in turbulence. Physical Review Fluids 7 (8), 083602\.
  * Rodriguez et al. (2011) Rodriguez, I., Borell, R., Lehmkuhl, O., Segarra, C. D. P. & Oliva, A. 2011 Direct numerical simulation of the flow over a sphere at re = 3700. Journal of Fluid Mechanics 679, 263–287.
  * Savva & Bush (2009) Savva, N. & Bush, J. W. M. 2009 Viscous sheet retraction. Journal of Fluid Mechanics 626, 211–240.
  * Shinjo (2018) Shinjo, J. 2018 Recent advances in computational modeling of primary atomization of liquid fuel sprays. Energies 11 (11), 2971.
  * Tang et al. (2021) Tang, K., Mostert, W., Fuster, D. & Deike, L. 2021 Effects of surface tension on the richtmyer-meshkov instability in fully compressible and inviscid fluids. Physical Review Fluids 6 (11), 113901.
  * Theofanous (2011) Theofanous, T. G. 2011 Aerobreakup of newtonian and viscoelastic liquids. Annual Review of Fluid Mechanics 43, 661–690.
  * Theofanous et al. (2012) Theofanous, T. G., Mitkin, V. V., Ng, C. L., Chang, C. H., Deng, X. & Sushchikh, S. 2012 The physics of aerobreakup. ii. viscous liquids. Physics of Fluids 24 (2), 022104\.
  * Troitskaya et al. (2017) Troitskaya, Y., Kandaurov, A., Ermakova, O., Kozlov, D., Sergeev, D. & Zilitinkevich, S. 2017 Bag-breakup fragmentation as the dominant mechanism of sea-spray production in high winds. Scientific Reports 7 (1), 1–4.
  * Troitskaya et al. (2018) Troitskaya, Y., Kandaurov, A., Ermakova, O., Kozlov, D., Sergeev, D. & Zilitinkevich, S. 2018 The ”bag breakup” spume droplet generation mechanism at high winds. part i: Spray generation function. Journal of Physical Oceanography 48, 2168–2188.
  * Vassallo & Ashgriz (1991) Vassallo, P. & Ashgriz, N. 1991 Satellite formation and merging in liquid jet breakup. Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 433 (1888), 269–286.
  * Veron (2015) Veron, F. 2015 Ocean spray. Annual Review of Fluid Mechanics 47, 507–538.
  * Villermaux & Bossa (2009) Villermaux, E. & Bossa, B. 2009 Single-drop fragmentation determines size distribution of raindrops. Nature Physics 5 (9), 697–702.
  * Vledouts et al. (2016) Vledouts, A., Quinard, J., Vandenberghe, N. & Villermaux, E. 2016 Explosive fragmentation of liquid shells. Journal of Fluid Mechanics 788, 246–273.
  * Wang et al. (2018) Wang, Y., Dandekar, R., Bustos, N., Poulain, S. & Bourouiba, L. 2018 Universal rim thickness in unsteady sheet fragmentation. Physical Review Letters 120 (20), 204503.
  * Yang et al. (2017) Yang, W., Jia, M., Che, Z., Sun, K. & Wang, T. 2017 Transitions of deformation to bag breakup and bag to bag-stamen breakup for droplets subjected to a continuous gas flow. International Journal of Heat and Mass Transfer 111, 884–894.
  * Yang et al. (2016) Yang, W., Jia, M., Sun, K. & Wang, T. 2016 Influence of density ratio on the secondary atomization of liquid droplets under highly unstable conditions. Fuel 174, 25–35.
  * Young (1995) Young, V. 1995 Liquid rocket engine combustion instability, , vol. 169. AIAA.
  * Zhang et al. (2019) Zhang, J., Chen, L. & Ni, M.-J. 2019 Vortex interactions between a pair of bubbles rising side by side in ordinary viscous liquids. Physical Review Fluids 4 (4), 043604.
  * Zhao et al. (2011) Zhao, H., Liu, H., Xu, J. & Li, W. 2011 Experimental study of drop size distribution in the bag breakup regime. Industrial and Engineering Chemistry Research 50, 9767–9773.
  * Zhou et al. (2021) Zhou, Y., Williams, R. J. R, Ramaprabhu, P., Groom, M., Thornber, B., Hillier, A., Mostert, W., Rollin, B., Balachandar, S., Powell, P. D. & others 2021 Rayleigh–taylor and richtmyer-meshkov instabilities: A journey through scales. Physica D: Nonlinear Phenomena p. 132838.
  * Zotova et al. (2019) Zotova, A. N., Troitskaya, Y. I., Sergeev, D. A. & Kandaurov, A. A. 2019 Direct numerical simulation of bag-breakup-mechanism of sea spray generation in strong winds. Journal of Physics: Conference Series 1163 (1), 012028.