Bag Film Breakup of Droplets in Uniform Airflows

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).

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.

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,

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,

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,

As $\Delta$ is the smallest length scale at which necks of thinning filaments can be represented, $L$ sets the length scale at which liquid filament breakup occurs.

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$, which defines the critical thickness $h_{c}\equiv 3D/2^{L_{\rm sig}}$, the smallest length scale at which liquid films can be presented as below it they will be artificially perforated by the MD algorithm. The signs of the computed quadratic moments indicate the local shape of the interface. If a film with thickness not larger than $h_{c}$ 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 with 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 requires the grid resolution level satisfy $L\geq 12$ for simulations with domain size $D=15d_{0}$. The highest grid resolution level we set in our present 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 diameter $d_{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.

We first discuss the initiation period of aerobreakup, defined by Jackiw & Ashgriz (2021) as $0\leq t\leq T_{i}$, where $T_{i}$ is the time when the droplet reaches its minimal streamwise thickness. 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.

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 can 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 a result of RT Instability developing on the windward surface due to wind acceleration. Li et al. (2019) predicted a critical instantaneous Bond number $Bo_{c}\equiv\rho_{l}\alpha d_{0}^{2}/4\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 windward side of the periphery 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 more 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.

For validation of our numerical results, we present in Fig. 3 the evolution of the maximum spanwise radius of the drop $R_{m}$ and the streamwise length of the bag $L_{\rm bag}$ measured from our axisymmetric and 3D numerical simulations at $We=15,\,Oh=2.5\times 10^{-3}$, and compare them with the available experimental results of Jackiw & Ashgriz (2021) and Flock et al. (2012). It can first be seen that the axisymmetric and 3D numerical results agree excellently until $t\approx 1.5\tau$, when the axisymmetric simulation shows a smaller bag length in fig. 3b. This late-time deviation most likely arises from the lack of 3D flow instability development in axisymmetric simulations (Marcotte & Zaleski, 2019), which may break the symmetry of the bag and limits its streamwise growth. Both our axisymmetric and 3D simulation results agree well with the experimental data of Jackiw & Ashgriz (2021) up to $t/\tau=1$, after which the experimental results show faster growth in both $R_{m}$ and $L_{\rm bag}$. This may be due to the sensitivity of the flattened drop to difference in the ambient flow conditions, as in our numerical simulations the air-phase flow remains laminar, whereas the experimental configuration of Flock et al. (2012) and Jackiw & Ashgriz (2021) in fact produces air-phase turbulence, which has been shown by Zhao et al. (2019) to be capable of increasing the height and width of bags at late time (see e.g., their fig. 6). More specifically, Jackiw & Ashgriz (2021) used a 5-gauge air needle (whose diameter $D_{n}$ is only 2.48 times of the droplet diameter $d_{0}$) to generate air jets with centreline Reynolds number $Re_{a}$ of $5.2\times 10^{3}\sim 2.5\times 10^{4}$, apart from needles for suspending the drop within such air jets. In the case of Flock et al. (2012), air jets are produced through a nozzle with diameter $D_{n}\approx 11d_{0}$, but the airflow is also turbulent with $Re_{a}=1.8\times 10^{4}$. The results of Jackiw & Ashgriz (2022) are obtained from single experimental runs without being ensemble-averaged, which may lead to larger variations in their results, as also noted in the comparison of numerical results by Ling & Mahmood (2023). Additionally, note that in Fig. 3a, the experimental results of Jackiw & Ashgriz (2021) and Flock et al. (2012) show some mutual disagreement in the spanwise radius values within the range of $t/\tau\leq 1.5$.

The bag lengths and widths recorded at various $Oh$ values at the point of breakup are also included as scattered points in fig. 3, which we will return to in §4.4. It can be seen that our bags approach breakup within the time range of $1.74\leq t/\tau\leq 1.91$, earlier than the experimental results of Jackiw & Ashgriz (2021) ($t=2.2\tau$ for $We=15.3$). On the other hand, Flock et al. (2012) did not report the exact time at which bag breakup is initiated. This earlier breakup time is associated with the limit of grid resolution, and hence $L_{\rm sig}$ on the critical thickness at which the bag film is perforated by the MD algorithm, as at $L_{\rm sig}=13$, the critical thickness is $3D/2^{L_{\rm sig}}=5.5\times 10^{-3}d_{0}$, which is a few times larger than the experimental value of $h/d_{0}=1.2\times 10^{-3}$ as found by Jackiw & Ashgriz (2022), and $5\times 10^{-5}\leq h/d_{0}\leq 5\times 10^{-4}$ by Opfer et al. (2014). This is a limitation present in all numerical simulations of droplet aerobreakup, as is also noted in the recent work of Ling & Mahmood (2023).

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

where $a$ is the axial stretching rate near the frontal stagnation point and approximated as $a\simeq 6$; $\tau$ is the characteristic deformation time introduced in §1; and $T_{\rm bal}$ is the time when a constant streamwise deformation rate is reached, taken as $0.125\tau$ according to the experimental results (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.

We now investigate (6) using our axisymmetric numerical results. Figs 4(a) and (b) show respectively the influence of $We$ and $Oh$ on the measured instantaneous spanwise growth rate $\tilde{\dot{R}}$, where the tilde indicates normalisation by the theoretical value (6). We also include the growth rate evolution obtained by numerically differentiating the experimental data presented in Fig. 28 of Jackiw & Ashgriz (2021) for comparison. For the small $Oh$ value of $10^{-3}$, fig. 4a indicates that the spanwise growth rate $\tilde{\dot{R}}$ reaches a plateau with relatively small variations around $t=0.3\tau$, where the prediction of (6) matches qualitatively with the measured $\tilde{\dot{R}}$ values. We note that while Jackiw & Ashgriz (2021) set $T_{\rm bal}=0.125\tau$ as an a posteriori estimation based on the evolution of $R_{m}$ rather than $\tilde{\dot{R}}_{m}$ when analysing their Fig. 17(b), our results agree well with the spanwise growth rate computed from their experimental data up to $t=0.74\tau$, with their data also reaching a plateau around $t=0.3\tau$. The growth rate of Jackiw & Ashgriz (2021) becomes much larger than ours for $t>0.74\tau$, corresponding to the larger $R_{m}$ values observed in fig. 3a, which is possibly a result of air-phase turbulence as previously discussed. For cases at $Oh=10^{-3}$, this period of constant $\tilde{\dot{R}}$ ends around $t=0.55\tau$, after which $\tilde{\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. 4b suggests that as $Oh$ increases beyond $2.5\times 10^{-3}$, the late-time peaking of $\tilde{\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 $\tilde{\dot{R}}$ is observed, the liquid is being pushed from the frontal surface to the windward side of the periphery, where the maximum spanwise radius is reached. However, at $t=0.6\tau$ when the peaking behaviour is observed, a bulge appears downstream and causes a location shift where the maximum spanwise radius $R$ is reached. This bulging behaviour is also present in the growth rate evolution computed from the experimental data of Jackiw & Ashgriz (2021) at $Oh=2.7\times 10^{-3}$, but virtually absent when $Oh=10^{-2}$ in our numerical simulations, as shown in fig. 2b, where the periphery of the droplet contour only flattens over time.

To provide insights into physical mechanisms governing the peak in the spanwise growth rate observed for low Oh values in fig. 4, we plot in fig. 5 the pressure distribution and streamlines near the drop periphery, when the peaks in the spanwise growth rate $\tilde{\dot{R}}$ are reached in fig. 4b 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 windward side of the periphery, 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 $\tilde{\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. 5a, 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 $\tilde{\dot{R}}$ after the peak. Notably, the droplet contour in fig. 5b at $Oh=10^{-2}$ lacks craters at the axis of symmetry and bulges at the periphery, and therefore more closely resembles the cylindrical shape of the deforming drop assumed in the derivations of (Jackiw & Ashgriz, 2021), which may explain why the match with the inviscid model (6) is improved as $Oh$ is increased.

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 6a shows the air- and liquid-phase pressure on either side of the drop surface as functions of the arclength $l$ traversing along the axisymmetric droplet contour in the clockwise direction. 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/d_{0}\leq 0.6$, which corresponds to the windward face of the drop; whereas the profile at $l/d_{0}>0.6$ deviates from the potential-flow solution due to flow separation, characterised by a second minimum around $l/d_{0}=1.2$. 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 $4\sigma/d_{0}$. As the droplet flattens over time, the air-phase pressure profile on the windward surface increases and the first minimum moves upstream, deviating from the potential-flow solution. In the meantime, the change in liquid-phase pressure for $l/d_{0}\leq 0.37$ is relatively small, and the air- and liquid-phase pressure profiles cross over each other at $l/d_{0}=0.37$ 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/d_{0}=0.85$ observed at $t/\tau=0.4$ becomes a maximum at $t/\tau=0.6$, which corresponds to the bulge formation observed in fig. 5a that leads to the deviation from (6) (Jackiw & Ashgriz, 2021).

Figure 6b shows the airphase axial velocity $u_{x}$ 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 $4/\pi$, the extreme values corresponding to spherical and pancake drop shapes as noted in Jackiw & Ashgriz (2021); crossing over the intermediate value of $2\sqrt{2}$ proposed by Kulkarni & Sojka (2014). The decrease of $a$ corresponds to the air-phase pressure increase on the windward surface as observed in fig. 6a via the following equation:

Consequently, we conclude that our numerical results reproduce the prediction (6) of Jackiw & Ashgriz (2021) that there exists a period characterised by a constant spanwise radius growth rate $\tilde{\dot{R}}$. Furthermore, we find that the later deviation from (6) is characterised by a peak in $\tilde{\dot{R}}$, which is caused by the capillary deceleration of liquid influx into the drop periphery that causes bulge formation at low $Oh$ values. The increase in $Oh$ eliminates the bulge and the frontal crater on the droplet surface, and the droplet acquires 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.

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. 7a 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 7a 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 7b, 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.

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$,

where the exponential decay rate $\lambda$ is given as $4/\tau$. We compare our numerical results with the predictions of model (8) in fig. 8. The film 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. 8, the exponential decay rate $\lambda$ is initially 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. 7. 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 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.

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

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 number of small fragments without a well-defined fragmentation mechanism; fig. 9a 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 the parameters of the MD algorithm can be adjusted so that sufficiently few holes are formed in these simulations and this mass loss becomes insignificant (see below).

The holes created by the MD perforation mechanism resemble those appearing in experiments (such as Lhuissier & Villermaux (2013); Jackiw & Ashgriz (2022); Kant et al. (2022)) - see fig. 9b. 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/or 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 first demonstrate the effects of the MD algorithm on the breakup behaviour of the bag film in fig. 9. The film rupture behaviour is qualitatively different with and without application of the MD algorithm. Figure 9a 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. 9b and 9c 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 these may collide with one another and form a few long stretching liquid bridges (Agbaglah, 2021), which are distinct from the numerous short and irregular bridges observed with VOF breakup. Note that fig. 9b still shows some VOF breakup behaviour which is absent in fig. 9c. 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. 9b, $\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. 9c.

Figure 10 further compares the grid convergence behaviour for the size and velocity distribution of bag fragments without and with the MD algorithm applied. The fragment data are sampled at different times throughout the bag film breakup period, and then collected and binned based on the equivalent fragment diameter $d$ to produce the size and velocity distribution functions. While the distribution functions presented in figs. 10a-10d are for single realisations, those in figs. 10e-10f are ensemble-averaged for each bin over different realisations; and we have verified that the total number of bins does not significantly influence the shape of size and velocity distributions. Figure 10a 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. 10c 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. 10a 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 by Chirco et al. (2022); and size distributions at different grid and signature levels appear to overlap for $d\geq 8\Delta_{\rm sig}$ despite these scatters. Figure 10e further presents the ensemble-averaged size distribution functions obtained with the MD algorithm applied, which features much smaller range of scatter, as indicated by the confidence bounds represented by the grey shade at $L=14,\,L_{\rm sig}=13$, showing clearly that the distributions of fragment statistics overlap for $d\geq 8\Delta_{\rm sig}$ at different values of $L$ and $L_{\rm sig}$. From this we conclude that the ensemble-averaged data are grid-converged for $d\geq 8\Delta_{\rm sig}$ and $13\leq L\leq 14,L>L_{\rm sig}$. Moreover, together figs. 10a,  10c and  10e establish that the lack of grid convergence in the no-MD case (fig. 10a) is attributable not to the scatter of individual realisations, but specifically to numerically uncontrolled VOF-breakup.

We also include the experimental data of Guildenbecher et al. (2017) obtained with $We=13.8$ and $Oh=5.43\times 10^{-3}$ in fig. 11, together with exponential and log-normal models fitted to their data. Guildenbecher et al. (2017) noted the difference between the size distributions obtained using two experimental techniques with different resolution levels, and expressed most confidence in the upper tail of the distributions satisfying $d\geq 0.01d_{0}$. Within this size range, the tails of our size distribution and those of Guildenbecher et al. (2017) show excellent agreement, 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. Both the exponential and the log-normal model are found to match well with the current size distribution functions for $d\geq 8\Delta$, while both differ from the current results within the range of $d<8\Delta$, which may suggest that no single function can represent the complete spectrum of the current size distribution of bag fragments.

It is noted that in fig. 10e, the fragment statistics are not fully converged for $d\leq 8\Delta_{\rm sig}$, where compared with its counterparts at $L=13$, the size distribution at $L=14$ shows more fragments satisfying $\Delta\leq d\leq\Delta_{\rm sig}$, and fewer fragments with $\Delta_{\rm sig}\leq 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). On the other hand, the size distribution function of very large fragments satisfying $d\geq 0.1d_{0}$ show relatively larger range of scatter compared with their smaller counterparts around $8\Delta_{\rm sig}$, which likely arises from the smaller number of these fragments produced in each ensemble realisation and can be further reduced by increasing the ensemble size. 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 10b, 10d and 10f 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. 10f, which is not observed in fig. 10b where large scatters across various grid levels are present. Interestingly, in fig. 10f, 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$.

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 (Vledouts et al., 2016). As is noted in §3.1, our film is thicker and breaks up earlier compared with experimental results due to the limit of grid resolution. This leads 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 results show that many interesting breakup mechanisms can already be captured with this choice of $L_{\rm sig}$, which we will present further below.

We first show in fig. 12 the time evolution of the instantaneous distributions of the size, axial and radial speed distributions of the fragments produced from bag breakup. Figure 12a 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$ that is very close to the ensemble- and time-averaged size distribution function. 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 and the fragment size distribution is well-represented by time-averaged results. The remaining rim will then disintegrate at still later times, whose investigation we leave for future work.

Figures 12b and 12c 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, with the ensemble-wide variations of velocity components in each averaging bin shown in error bars. 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/d_{0}=5\times 10^{-3}$ in the distributions of both $u_{x}$ and $u_{r}$, whose location does not appear to vary with time. Moreover, despite the presence of velocity variations, fig. 12b 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. 12c. 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 velocity component of the 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.

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. 13 shows the breakup of a stretched long ligament neighbouring two enlarging holes into a series of small drops. As the ligament is itself connected to the main drop, there is a significant size difference between the parent and child drops produced from its breakup, which is an example of non-local breakup events (see Eq. (11) in §4.3 for a definition of non-local breakup). It can be seen from fig. 13a 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. 13b, and forms a primary drop which continues to undergo periodic prolate-oblate shape oscillations resembling droplets produced by breaking Rayleigh jets (Hu et al., 2021), as highlighted in the red boxes in figs. 13b-13h. 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 debris, and then split into large primary and small satellite drops via the well-known end-pinching mechanism (Castrejón-Pita et al., 2012; Pal et al., 2021), which is an example of local breakup events as the parent (the elongated debris) and child (satellite drops) do not differ significantly in their sizes. Under certain circumstances, the primary and satellite drops 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 14 first shows an example of short ligament breakup and its eventual contraction into a single droplet, as highlighted in the blue boxes. Compared with the breakup of long ligaments demonstrated in fig. 13, 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.

Another type of fragmentation mechanism can also be identified in fig. 14; 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. 12a; namely, that larger droplets are relatively few, and produced at generally later times during the fragmentation process.

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.14, 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}=1.2d_{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 $2d_{0}$, as shown in both our numerical results and the experimental data of Jackiw & Ashgriz (2021) (see e.g. their fig. 30). Figure 15 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.

Lastly, we observed a few examples showing the destabilisation of receding liquid rims, as demonstrated in fig. 16. 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 the type of hydrodynamic instability involved here; we observe in fig. 16d 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. 16). 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 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.

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 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).

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 highly 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 17a 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. 17a as a function of $d_{p}$, which is defined as follows:

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.05d_{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.05d_{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.05d_{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.05d_{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. 17b 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. 17a, 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}$:

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. 17b 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. 17b, $\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. 17b, we introduce the following criteria for determining ‘non-local’ breakup and coalescing events, respectively:

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. 17b, 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).

We further plot in fig. 18 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 exist ‘multistep’ breakup events in our aerobreakup simulation outputs, our results in fig. 18 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.

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. 19a shows the evolution of surface energy $E_{s}$ of individual small fragments recorded from simulations run at $L=14$ and $L_{\rm sig}=13$, with the records of only a few representative fragments highlighted for clarity. The steady-state values of surface energy is also computed based on the volume of corresponding fragments, and plotted in grey dashed lines for reference. It is seen that the oscillation frequency and amplitude vary for each fragment, but all of them clearly demonstrate decaying oscillation behaviour, 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 these small fragments at two different grid resolution configurations ($L=14,\,L_{\rm sig}=13$ and $L=13,\,L_{\rm sig}=11$), and plot them against the fragment radius in fig. 19b, where an excellent agreement is found for small fragments within the diameter range of $0.01d_{0}\leq d\leq 0.1d_{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:

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. 13. Both the viscous and inviscid theoretical model of Prosperetti (1980) are plotted in fig. 19, 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. 19b, 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 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).

Having provided an overview of the physical mechanisms governing bag breakup and the subsequent evolution patterns of the fragments at a specific low $Oh$ value of $10^{-3}$ in §4.2, in this section we increase $Oh$ up to 0.05, 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).

We first provide in fig. 20 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 is first seen that as $Oh$ increases from very low ($Oh\leq 0.005$) to moderate ($0.005\leq Oh\leq 0.05$) values, 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 moderate-$Oh$ 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. 16, 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 their formation out of colliding hole rims, and then break up into ‘primary’ and ‘satellite’ drops that differ significantly in their sizes. When $Oh$ increases to the moderate value of $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 their free ends 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, 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.

Figure 21 further shows the evolution of the number of fragments satisfying $d\geq\Delta_{\rm sig}$, the dependence of bag length and width at the onset of bag breakup on the $Oh$ values, and the time- and ensemble-averaged size and speed distribution functions for the three specific $Oh$ values selected in fig. 20. Figure 21a 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 feeds 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}$. The length $L_{\rm bag}$ and width $d_{\rm bag}$ of the bags just before their breakup for different $Oh$ values are measured and shown in fig. 21b (which were also presented in fig. 3 as scattered points for validation of our numerical results), where it can be seen that our bags are ‘flattened’ in shape (satisfying $d_{\rm bag}>L_{\rm bag}$), and that for $Oh\geq 5\times 10^{-4}$, both $L_{\rm bag}$ and $d_{\rm bag}$ decrease with increasing $Oh$, suggesting that the bag area indeed becomes smaller as $Oh$ increases.

Figure 21c 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. 20l) 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. 21d 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 ambient airflow becoming much less turbulent as its viscosity also increases under a fixed viscosity ratio $\mu^{*}$.

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. 22. It should be noted that due to the significant runtime required on supercomputers, the results at a few $Oh$ values in fig. 22 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. 22a, it is seen that $\bar{t}_{b}$ first remains almost independent of $Oh$ as the latter increases from the low value of $10^{-4}$ to the moderate value of $10^{-2}$, which is likely because the wake flow remains separated from the drop surface, and the thinning process of the bag before the onset of its breakup is determined by capillary and inertial effects, as discussed in §3.2. When $Oh$ exceeds the moderate value of $10^{-2}$, $\bar{t}_{b}$ is found to increase exponentially with $Oh$ as shown by the fitted model, which may be the consequence of both the transition of the wake region from a 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$:

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 film 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. 22b, $\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 bag films continues thinning 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, non-dimensionalised by $h_{f}$, match the order of magnitude of the thin-film fragment statistics studied by Kant et al. (2022), which are also included in fig. 22b 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$, which corresponds to the drop $Oh$ value exceeding the moderate value of 0.01. 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. 20 and fig. 21, 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.