<figure><img src="image_1.png" /></figure>

HAL
open science

# Drop Deformation Dynamics and Fragment Characteristics in the Bag-Breakup Regime

Kaitao Tang, Thomas Adcock, Wouter Mostert

► To cite this version:

Kaitao Tang, Thomas Adcock, Wouter Mostert. Drop Deformation Dynamics and Fragment Characteristics in the Bag-Breakup Regime. the 23rd Australasian Fluid Mechanics Conference, Dec 2022, Sydney, Australia. hal-03893409

HAL Id: hal-03893409

https://hal.science/hal-03893409

Submitted on 11 Dec 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Drop Deformation Dynamics and Fragment Characteristics in the Bag-Breakup Regime

Kaitao Tang¹*, Thomas Adcock¹ and Wouter Mostert¹

¹Department of Engineering Science, University of Oxford, Oxford OX1 3PJ, UK

*Email: kaitao.tang@eng.ox.ac.uk

## Abstract

We investigate numerically the aerobreakup of liquid droplets with a focus on the bag breakup regime. We solve the two-phase incompressible Navier-Stokes equation using the adaptive mesh refinement (AMR) technique with volume-of-fluid (VOF) interface reconstruction. To minimise the influence of mesh-induced breakup on fragment statistics which has prevented previous studies from reaching grid convergence, we utilise a recently developed numerical algorithm which artificially perforates thin fluid structures once their thickness decreases to a critical value independent of the grid size. We first show good agreement with theoretical predictions and clarify the influence of the surrounding airflow on the drop deformation. Next, we show grid convergence for fragment statistics and identify the physical mechanisms governing bag breakup.

## 1 Introduction

Liquid atomisation refers to the process where a bulk volume of liquid deforms and disintegrates into sprays, typically featuring various shapes and spanning large ranges of size and velocity [1]. Liquid atomisation is involved in a wide range of physical processes, including ocean-atmosphere interactions [2], precipitation and rain-drop dynamics [2, 3], pharmaceutical spray generation and pathogen transmission [4], etc.; and it has recently been shown that the fragmentation of bag-shaped surface ripples dominates sea spume generation under extreme wind conditions, producing large droplets with typical sizes of $10^2 \sim 10^3 \mu m$ [5] where the currently-available sea-spray generation functions (SS-GFs) show large range of scatter [2]. However, since the physics governing the fragmentation of bag films have not yet been firmly established, their influence on SSGFs have been difficult to quantify. Improving this understanding is the primary motivation of the present work.

As one of the simplest examples of liquid atomisation, the aerobreakup of liquid drops is typically modelled by placing a droplet with density $\rho_l$, viscosity $\mu_l$, diameter $d_0$ and surface tension $\sigma$ in an ambient gas flow with density $\rho_a$, viscosity $\mu_a$ and uniform incoming velocity $U_0$. Four non-dimensional controlling parameters can be proposed [3]:

$$ 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}. \qquad (1) $$

Among these, the Weber (We) and Ohnesorge number (Oh) respectively quantify the ratios of inertial to capillary forces, and viscous to capillary forces; and $\rho^*$ and $\mu^*$ are respectively the density and viscosity ratios of the liquid and gas phase. For sufficiently small We, the droplet oscillates without breakup, up to a critical value $We_c = 11 \pm 2$ [6] which marks the transition to the bag breakup regime, and is independent of viscous effects for $Oh \le 0.1$.

In the bag-breakup regime, the droplet first flattens and evolves into an oblate disc. Afterwards, a bag is blown out from the centre, thinning and approaching breakup, while the disc periphery forms a toroidal rim. There are several models describing the early deformation of the droplet [1, 6], but

<figure><img src="image_2.png" /></figure>
the internal flow mechanism first proposed by Villermaux and Bossa [3] and significantly developed by Jackiw and Ashgriz [7, 8] among others has been very promising. In this model, the droplet is assumed to maintain a spheroidal shape governed by axisymmetric and inviscid internal flow, with a well-defined air-phase stagnation point [1, 3, 7]. Recently, Jackiw and Ashgriz [7] used a simple semi-analytical argument to propose a constant growth rate of the drop spanwise diamater $R_m$, which has shown good agreement with experimental results [7]. However, this approach is somewhat empirical and cannot account for the complex interaction between wake vortices and drop surface [9]. The late-stage 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 [1]; and it has only recently been clarified [8] that the major pathways leading to bag fragmentation are the destabilisation and collision of hole rims as they recede over the curved bag and experiences centripetal acceleration. The size and velocity distribution data of aerobreakup fragments are still scarce [10], and it remains unclear how each breakup mechanism contributes to the fragment statistics given the large span in time and length scales.

While most of the early aerobreakup studies are experimental, recent development of computa-
tional power has enabled numerical investigations of aerobreakup; however, they are highly challeng-
ing due to the multiscale nature of the problem. 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
[11]. As a result, numerical convergence with respect to bag fragment statistics has not previously
been obtained to our knowledge. Using a controlled and well-defined perforation scheme to obtain
numerical convergence is therefore one of the major goals of this study.

In this work, we investigate both the early-time deformation and the late time fragmentation pro-
cess of a droplet in bag breakup. We first introduce in §2 the problem configuration and the numerical
method. We then compare our axisymmetric simulation results with previous theoretical predictions
in §3 for the early-time deformation period, and discuss in §4 the 3D grid convergence study of
fragment statistics and the physical mechanisms governing bag breakup.

## 2 Formulation and Methodology

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sketch showing the initial configuration of axisymmetric (a) and three-dimensional (b) droplet aerobreakup simulations. For the axisymmetric simulations, the axis of symmetry is located at the bottom.</figcaption></figure>

We consider both axisymmetric and three-dimensional (3D) problem geometries, and the corre-
sponding flow configurations are shown in Figs. 1a and 1b, respectively. To eliminate the influence
of finite domain size on aerobreakup behaviour, the domain width $D$ is set as $15d_0$. A zero-gradient
velocity boundary condition and a uniform incoming velocity $U_0$ are imposed respectively on the right
and left boundaries, while reflecting conditions are applied at the other boundaries. This initialisation
results in an impulsive droplet acceleration at the first time step, and induces a flow field satisfying both the incompressibility constraint and the conservation of linear momentum [9]. We fix $\rho^*$ and $\mu^*$ as 833 and 55 which are typical of air-water systems. We and Oh are varied within $12 \le We \le 20$ and $0.001 \le Oh \le 0.01$, respectively, which covers the bag and bag-stamen breakup regimes.

We use the open-source software Basilisk [12] to solve the following Navier-Stokes equations for two-phase incompressible, immiscible and isothermal flows,

$$ \nabla \cdot \mathbf{u} = 0, \qquad (2) $$

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \nabla \cdot \left[ \mu (\nabla \mathbf{u} + \nabla \mathbf{u}^T) \right] + \sigma \kappa \delta_s \mathbf{n}. \qquad (3) $$

Equations (2) and (3) are respectively the continuity and momentum equation, where $p$ is the fluid pressure. $\sigma\kappa\delta_s\mathbf{n}$ incorporates surface tension effects, where $\kappa$ and $\mathbf{n}$ are respectively the local interfacial curvature and normal vector, and the Dirac delta $\delta_s$ is non-zero only on the interface [13]. Compressible and mixed compressible-incompressible problems can also be considered using the appropriate Basilisk solver [14], but this is out of the scope of the present work.

The geometric volume-of-fluid (VOF) interface reconstruction method is applied, which solves the following advective equation,

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

where the VOF function $f$ equals 1 and 0 in the liquid and gas phase, respectively. $\delta_s\mathbf{n}$ in Eq. (3) is approximated as $\nabla f$ [15, 16], and $\kappa$ is calculated using the appropriate numerical derivatives of the interface geometry, which is reconstructed through the use of height functions [16]. The quad/octree adaptive mesh refining (AMR) scheme is adopted to reduce the computational cost at high resolution levels $L$ [17], which is defined using the minimum grid size $\Delta = D/2^L$. Our axisymmetric simulations are conducted on $L = 14$, while the 3D simulations are conducted on $L = 13$ and 14.

To achieve grid convergence of fragment statistics, we adopt and parallelise the manifold death (MD) algorithm of Chirco et al. [11], which randomly perforates thin films once their thickness decreases to a prescribed value independent of $\Delta$. This is realised by periodically computing quadratic moments of $f$ indicating the local interfacial topology at a given signature level $L_{\text{sig}} \le L$. When the detected local thickness reduces to $h_c = 3D/2^{L_{\text{sig}}} = 3\Delta \cdot 2^{L-L_{\text{sig}}}$, the algorithm randomly creates cavities on the film by modifying $f$ with a prescribed probability $p_{\text{perf}}$. While this approach changes the total fluid mass, the MD algorithm minimises this side effect by creating cavities with minimum sizes that allow for expansion, and limiting the maximum number of holes $n_H$ perforated within a calling interval $\Delta t_c$.

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 \equiv \sqrt{\rho_l/\rho_a d_0/U_0}$ [18] 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 results in the remainder of this study unless otherwise specified.

## 3 Early-time deformation

To provide an overview of the early-time droplet deformation process, we first present in Fig. 2a the development of droplet contours for $We = 15$, $Oh = 0.001$. It is found that the windward surface of the droplet continues moving downstream and pushing the liquid to the drop periphery, leading to gradual spanwise flattening; in the meantime a dimple develops on the windward surface that moves radially inwards and eventually evolves into a crater on the axis of symmetry. The leeward side of the droplet remains relatively stationary for a long time after an initial recession to the upstream.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a): Early-time development of droplet contours for axisymmetric simulation case with We = 15. Oh = 0.001; (b): Gas-phase pressure profiles around the droplet from axisymmetric and 3D cases with We = 15, Oh = 0.025. Dashed arrows indicate the profile evolution as time elapses.</figcaption></figure>

Figure 2b shows the gas-phase pressure as a function of the arc length *l* along the drop surface, where We = 15 and Oh = 0.025, which matches the rigid sphere case in [19] with a Reynolds number of 300. At very early time the pressure profile agrees in trend with the steady-state result of [19], especially in the wake region; and the axisymmetric and 3D results match very well over the entire deformation period. As the droplet flattens over time, the pressure profiles deviate from the rigid-sphere case, which is most apparent at the windward face of the drop.

Recently, Jackiw *et al.* [7] found experimentally a period of constant droplet spanwise growth rate *Ṙ*, and derived the following model assuming axisymmetric and inviscid flow,

$$ \dot{R} = \frac{R_0 T_{bal}}{4\tau^2} \left( a^2 - \frac{128}{We} \right), \qquad (5) $$

where the axial stretching rate *a* is approximated as 6, following the potential flow solution for a sphere; and *T*<sub>bal</sub> is the time when a constant streamwise deformation rate is reached, taken as 0.125τ according to the experimental results of [7].

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Development of instantaneous spanwise growth rate *Ṙ* at various We with Oh = 0.001 (a), and at various Oh with We = 15 (b). Numerical results are scaled using Eq. (5).</figcaption></figure>

Figures 3a and 3b show the influence of We and Oh on *Ṙ*, respectively. For Oh = 0.001, Figure 3a shows that *Ṙ* reaches a plateau around *t* = 0.3τ, where Eq. (5) matches qualitatively with the measured
$\dot{R}$ values. For $Oh = 0.001$, this period of constant $\dot{R}$ ends around $t = 0.55\tau$, after which $\dot{R}$ reaches a peak around $t = 0.6\tau$, indicating a deviation from Eq. (5) which becomes more conspicuous as We increases. On the other hand, Figure 3b suggests that as $Oh$ increases, the peaking of $\dot{R}$ attenuates, while the match with Jackiw et al.'s model (5) is improved, which is particularly interesting as Eq. (5) is based on inviscid flow assumptions. The experimental data in [7] are obtained for $Oh = 0.0027$ where the peaking behaviour is not apparent, which partially explains why it is not noted therein.

Figure 4 further shows that the gas flow separates from the droplet surface at the upper-left tip, creating a recirculating region with low pressure [9]. The streamlines within the droplet shows the liquid being pushed from the windward surface to the periphery driven by the air-phase pressure difference; and when the peak in $\dot{R}$ is reached for $We = 20$, $Oh = 0.001$, a liquid bulge appears around $x = 0.5$ in the recirculating region and causes a location shift where the maximum spanwise radius $R$ is reached, which could also be seen in Fig. 2a at $t \ge 0.6\tau$ for $We = 15$. Surface tension then induces a pressure increase in the bulge and decelerates the incoming liquid flow, thus causing the later decrease in $\dot{R}$. Notably, the droplet at $Oh = 0.01$ shows no bulging behaviour and much weaker gas-phase recirculation near the drop periphery in Fig. 4c, which suggests that the bulge formation and deviation from (5) is closely associated with the interactions between the drop surface and the wake vortices, which is in turn determined by liquid and gas viscosity.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Flow fields near the droplet before (a) and after (b) the peaks in $\dot{R}$ are reached for $We = 20$, $Oh = 0.001$, and at the peak in $\dot{R}$ for $We = 15$, $Oh = 0.01$ (c). Streamlines are shown in black with superimposed arrows indicating flow direction, while the background colour shows the flow pressure magnitude.</figcaption></figure>

# 4 Breakup of Bag Films

While there is still debate on the physical mechanism responsible for the onset of liquid film breakup, for bag films under normal acceleration, it has been argued that RT modulations arise across the film, causing perforation when the perturbation amplitude becomes comparable to the film thickness $h$ [3, 7]. In our simulations the bag film is perforated artificially to minimise the influence of numerical breakup, as discussed in §2; and here we will focus on bag dynamics after its nucleation.

We first show grid convergence for ensemble-averaged fragment statistics in Fig. 5 when applying the MD algorithm, where we set the perforation probability, calling interval and signature level as $p_{perf} = 1/17500$, $\Delta t_c = 0.25d_0/U_0$ and $L_{sig} = 12$ (see §2 for their definition). The film rupture behaviour is qualitatively different with and without application of the MD algorithm. Figure 5a is a snapshot for a simulation case run without using the MD algorithm, featuring spurious numerical breakup characterised by small-scale irregular corrugations and ligament breaking on the bag. Figures. 5b and 5c show that the MD algorithm is able to create large expanding holes, and reduce the influence of numerical breakup on fragment statistics. More specifically, Fig. 5b still shows some VOF breakup behaviour which is absent in Fig. 5c. This is because larger $L$ means smaller grid size $\Delta$ at which VOF breakup occurs, which postpones its onset and allows more time for the expansion
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Grid convergence study at grid level *L* = 12 and 13 for *We* = 15, *Oh* = 0.001. (a)-(c): Simulation snapshots showing fragmenting bag films at *t*/*τ* = 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*<sub>sig</sub> = 12. (d): Evolution of the instantaneous fragment size distribution; (e): ensemble- and time-averaged fragment size distribution for 1.906 ≤ *t*/*τ* ≤ 2.183.</figcaption></figure>

of holes created by the MD algorithm with $L_{\text{sig}}$ unchanged.

Figure 5d shows the evolution of the fragment size distribution for one ensemble realisation at grid
level *L* = 12 and 13, respectively; where each horizontal 'line' corresponds with a fragment with a
given size, formed at the time when the line first appears; variations in the radii of fragments corre-
spond roughly with their breakup or coalescence, although there is some numerical noise particularly
at small sizes. Initially only a few small fragments with radii close to the grid size Δ are produced,
followed by well-resolved larger fragments with radii *r* ≥ 4Δ where grid convergence is reached [20].
The small fragments close to grid size are mostly formed due to the breakup of stretched liquid lamella
bordering colliding holes [21], while larger fragments with *r* ≥ 4Δ are formed as 'nodes' adjacent to
three or more expanding holes detach from the bag [22], as shown in Figs 5b and 5c. Figure 5e
shows the ensemble- and time-averaged fragment size distribution function for *L* = 12 and 13, where
the statistics are approximately grid-converged for the tail with *r* ≥ 4Δ₁₂ when the MD algorithm
is used. In contrast, fewer fragments with *r* ≥ 4Δ₁₂ are generated without using the MD algorithm,
which is most likely because of the predominance of irregular tiny fragments characteristic of un-
controlled VOF breakup, as shown in Fig. 5a. It is noted that the difference between the *L* = 12 and
*L* = 13 data is most apparent for *r* ≤ 4Δ₁₃, where the *L* = 13 distribution extends further into the range
of *r* ≤ Δ₁₃. This is probably because the smallest fragments near the grid size are produced from cap-
illary breakup of slender ligaments, which are controlled by the grid resolution rather than the MD
algorithm. Lastly, the loss of liquid mass incurred by the MD algorithm does not exceed 0.023% for
*t*/τ ≤ 2.18 at both *L* = 12 and 13, showing minimal influence of MD on mass conservation.

Finally, we briefly discuss the types of fragmentation arising during the bag breakup. In Fig. 5,
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Snapshots showing destabilisation of the rim during hole expansion on the bag with We = 15, Oh = 0.001 (a-c). The direction of mean background airflow is from right to left.</figcaption></figure>

most fragments are formed by collision between receding rims, as in [21, 23], but other mechanisms
might also play a role. Figure 6 shows hole rim destabilisation for We = 15, Oh = 0.001, which to
our knowledge has not been shown in previous numerical studies. The simulation is run at L = L$_{sig}$ =
14. It could be seen that the rim bordering the major hole protrudes out from the bag surface as it
collects up liquid mass during receding, and regularly-spaced digitations develop on the upper half
of the rim as time elapses. Jackiw and Ashgriz [8] suggests that rims receding at the Taylor-Culick
velocity v$_{TC}$ ≡ $\sqrt{2σ/ρ_lh}$ over a bag with radius of curvature R are subject to the centrifugal inertial
force a$_{rr}$ = $v_{TC}^2/R$, which regulates the rim thickness b via the universal Bond-number principle [24]
Bo$_{rim}$ ≡ $ρ_l b^2 a/σ$ = 1. The experimental results in [8] further suggest that the Rayleigh-Plateau (RP)
instability governs rim digitation, while Lhuissier and Villermaux [25] proposes an RT instability
mechanism for bursting surface bubbles, which we aim to further investigate in our future work.

**5 Conclusions**

This study carries out axisymmetric and 3D numerical simulations of droplet aerobreakup, and
utilises a new film perforation algorithm to minimise pollution of fragment statistics by spurious
numerical breakup. Our results are validated by good agreement with [19]. We reproduce the con-
stant spanwise growth rate predicted in [7], and explain deviation from this model by the interaction
between the drop surface and the wake vortices. Grid convergence for fragment statistics and prelim-
inary results showing the destabilisation of receding rims are presented.

**Acknowledgements**

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**

[1] D. R. Guildenbecher, C. López-Rivera, and P. E. Sojka, Secondary atomization, *Exp. Fluids* **46**, 371–402 (2009).

[2] F. Veron, Ocean spray, *Annu. Rev. Fluid Mech.* **47**, 507–538 (2015).
[3] E. Villermaux and B. Bossa, Single-drop fragmentation determines size distribution of rain-drops, *Nat. Phys.* **5**, 697–702 (2009).

[4] L. Bourouiba, The fluid dynamics of disease transmission, *Annu. Rev. Fluid Mech.* **53**, 473–508 (2021).

[5] Y. Troitskaya, A. Kandaurov, O. Ermakova, D. Kozlov, D. Sergeev, and S. Zilitinkevich, The "bag breakup" spume droplet generation mechanism at high winds. part i: spray generation function, *J. Phys. Oceanogr.* **48**, 2168–2188 (2018).

[6] T. Theofanous, Aerobreakup of newtonian and viscoelastic liquids, *Annu. Rev. Fluid Mech.* **43**, 661–690 (2011).

[7] I. M. Jackiw and N. Ashgriz, On aerodynamic droplet breakup, *J. Fluid Mech.* **913** (2021).

[8] I. M. Jackiw and N. Ashgriz, Prediction of the droplet size distribution in aerodynamic droplet breakup, *J. Fluid Mech.* **940** (2022).

[9] F. Marcotte and S. Zaleski, Density contrast matters for drop fragmentation thresholds at low ohnesorge number, *Phys. Rev. Fluids* **4**, 103604 (2019).

[10] H. Zhao, H. F. Liu, J. L. Xu, and W. F. Li, Experimental study of drop size distribution in the bag breakup regime, *Ind. Eng. Chem. Res.* **50**, 9767–9773 (2011).

[11] L. Chirco, J. Maarek, S. Popinet, and S. Zaleski, Manifold death: a volume of fluid implementation of controlled topological changes in thin sheets by the signature method, *J. Comput. Phys.* **111468** (2022).

[12] S. Popinet, Basilisk flow solver and pde library, Available at: http://basilisk.fr, 2019.

[13] S. Popinet, Numerical models of surface tension, *Annu. Rev. Fluid Mech.* **50**, 49–75 (2018).

[14] K. Tang, W. Mostert, D. Fuster, and L. Deike, Effects of surface tension on the richtmyer-meshkov instability in fully compressible and inviscid fluids, *Phys. Rev. Fluids* **6**, 113901 (2021).

[15] J. U. Brackbill, D. B. Kothe, and C. Zemach, A continuum method for modeling surface tension, *J. Comput. Phys.* **100**, 335–354 (1992).

[16] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, *J. Comput. Phys.* **228**, 5838–5866 (2009).

[17] J. A. Van Hooft, S. Popinet, C. C. Van Heerwaarden, S. J. Van der Linden, S. R. de Roode, and B. J. Van de Wiel, Towards adaptive grids for atmospheric boundary-layer simulations, *Boundary-Layer Meteorol.* **167**, 421–443 (2018).

[18] J. Nicholls and A. Ranger, Aerodynamic shattering of liquid drops. *AIAA J.* **7**, 285–290 (1969).

[19] J. Magnaudet, M. Rivero, and J. Fabre, Accelerated flows past a rigid sphere or a spherical bubble. part 1. steady straining flow, *J. Fluid Mech.* **284**, 97–135 (1995).

[20] W. Mostert, S. Popinet, and L. Deike, High-resolution direct simulation of deep water breaking waves: transition to turbulence, bubbles and droplets production, *J. Fluid Mech.* **942** (2022).

[21] G. Agbaglah, Breakup of thin liquid sheets through hole-hole and hole-rim merging, *J. Fluid Mech.* **911** (2021).

[22] H. Lhuissier and E. Villermaux, ‘effervescent’ atomization in two dimensions, *J. Fluid Mech.* **714**, 361–392 (2013).

[23] B. Néel, H. Lhuissier, and E. Villermaux, ‘fines’ from the collision of liquid rims, *J. Fluid Mech.* **893** (2020).

[24] Y. Wang, R. Dandekar, N. Bustos, S. Poulain, and L. Bourouiba, Universal rim thickness in unsteady sheet fragmentation, *Phys. Rev. Lett.* **120**, 204503 (2018).

[25] H. Lhuissier and E. Villermaux, Bursting bubble aerosols, *J. Fluid Mech.* **696**, 5–44 (2012).