3.1. Rim node instability

3.1.1. The Rayleigh–Taylor instability

Zhao et al. (Reference Zhao, Liu, Li and Xu2010) suggest that the rim nodes are the result of a Rayleigh–Taylor instability on the rim, which has a maximum susceptible wavelength, $\lambda _{RT}$, of

(3.1)\begin{equation} \lambda_{RT} = 2 {\rm \pi}\sqrt{\frac{3 \sigma}{\rho_l a}}, \end{equation}

where $a$ is the acceleration of the deforming droplet. The acceleration $a$ can be estimated by considering the drag force on the flattened droplet as

(3.2)\begin{equation} a = \frac{3}{4} C_D \frac{U^2}{d_0} \frac{\rho_g}{\rho_l} \left( \frac{2R}{d_0} \right)^2, \end{equation}

where $C_D$ is the drag coefficient of the flattened droplet and $2R/d_0$ is the extent of the droplet deformation at the instant the instability takes hold. Upon substitution, the wavelength of the Rayleigh–Taylor instability on the liquid rim is

(3.3)\begin{equation} \frac{\lambda_{RT}}{d_0} = \frac{4 {\rm \pi}}{\sqrt{C_D We}} \left( \frac{d_0}{2R} \right). \end{equation}

Although this expression was derived to predict the instability wavelength over the whole droplet with the aim of modelling the Rayleigh–Taylor piercing breakup mechanism at very high $We$ (Theofanous et al. Reference Theofanous, Li and Dinh2004), Zhao et al. (Reference Zhao, Liu, Li and Xu2010) argued that the same analysis can also be used to predict the instability wavelength of the rim.

From (3.3), the two undetermined variables are $C_D$ and $2R$. Physically, the drag coefficient can vary between that of a sphere ($C_D=0.4$) and that of a disk ($C_D=1.2$). Since the highest acceleration will occur when the droplet takes the approximate shape of a disk, the value $C_D=1.2$ is practical for the case of the Rayleigh–Taylor instability.

Zhao et al. (Reference Zhao, Liu, Li and Xu2010) assume that the instability takes hold when the droplet reaches its maximum cross-sectional deformation, $R_{max}$, prior to bag blow-out, which is slightly after the rim is initiated. Zhao et al. (Reference Zhao, Liu, Li and Xu2010) provided the following correlation for $R_{max}$:

(3.4)\begin{equation} \frac{2R_{max}}{d_0} = \frac{2}{1 + \exp ({-}0.0019 We^{2.7})}. \end{equation}

In the analysis of Zhao et al. (Reference Zhao, Liu, Li and Xu2010), the correlation used for the prediction of $N$ is that of Hsiang & Faeth (Reference Hsiang and Faeth1992) ($2R_{max}/d_0 = 1 + 0.19 \sqrt {We}$); however, the correlation given by (3.4) provides a more consistent prediction of $N$, as is shown later in this section.

3.1.2. The Rayleigh–Plateau instability

The Rayleigh–Plateau instability wavelength, $\lambda _{RP}$, depends only on the thickness of the rim at its initiation, $h_i$ (i.e. the minor diameter of the toroidal rim), as (Rayleigh Reference Rayleigh1878)

(3.5)\begin{equation} \frac{\lambda_{RP}}{d_0} = 4.5 \frac{h_i}{d_0}, \end{equation}

where it is assumed that the rim is approximately a cylinder with periodic symmetry. The thickness $h_i$ is given by our previous model (Jackiw & Ashgriz Reference Jackiw and Ashgriz2021); however, in the present work, we consider a minor modification to the model, which is discussed in detail in Appendix A. The model for $h_i$ is given by

(3.6)\begin{equation} \frac{h_i}{d_0} = \frac{4}{We_{rim} + 5 \left( \dfrac{2R_i}{d_0} \right) - 4 \left( \dfrac{d_0}{2R_i} \right) } - 0.05, \end{equation}

where $We_{rim} = \rho _l \dot {R}^2 d_0 / \sigma$ is the rim $We$, for which $\dot {R}$ is given by

(3.7)\begin{equation} \frac{2\dot{R}}{d_{0}}=\frac{1.125}{\tau}\left(1-\frac{32}{9 We}\right). \end{equation}

The modification to the model is in the definition of the cross-sectional deformation at initiation, $2R_i$, which is now given by the following empirical relation (see (A1)):

(3.8)\begin{equation} \frac{2R_i}{d_0} = 1.63 - 2.88 \exp ({-}0.312 We). \end{equation}

Other combinations of the instability theories and extents of deformation are compared in Appendix B. The two presented here are the most reasonable and accurate of the possible combinations for both the Rayleigh–Taylor and Rayleigh–Plateau instabilities.

3.1.3. Instability breakup time

The first comparison of the plausibility of the Rayleigh–Taylor and Rayleigh–Plateau instabilities in the formation of the rim nodes is in the comparison of their characteristic breakup times, $t_b$. The breakup time of a liquid column is found by the time required for an initial perturbation to grow to the radius of the column, given by Grant & Middleman (Reference Grant and Middleman1966) and Ashgriz (Reference Ashgriz2011)

(3.9)\begin{equation} t_b = \frac{C_t}{\omega_{max}}, \end{equation}

where $C_t = \ln ( (h_i/2) / \zeta _0)$ is a constant related to the initial perturbation, $\zeta _0$, and the radius of the column, $(h_i/2)$, and $\omega _{max}$ is the growth rate of the most unstable wave. Since $\zeta _0$ is very small ($O({\rm nm})$), $C_t$ must be determined experimentally. As a point of reference, Grant & Middleman (Reference Grant and Middleman1966) give $C_t=13$ for the breakup of a laminar jet by the Rayleigh–Plateau instability.

The maximum susceptible growth rates for the Rayleigh–Taylor, $\omega _{max, RT}$, and Rayleigh–Plateau, $\omega _{max, RP}$, instabilities are given by

(3.10)\begin{equation} \omega_{max, RT} = \left( \frac{a^3 \rho_l}{2 \sigma} \right)^{1/4} \end{equation}

and

(3.11)\begin{equation} \omega_{max, RP} = 0.97 \sqrt{\frac{\sigma}{\rho_l h_i^2}}, \end{equation}

respectively (Chandrasekhar Reference Chandrasekhar1961). The acceleration, $a$, is estimated using (3.2). Additionally, the value of $2R_{max}$ used in the Rayleigh–Taylor theory is calculated using (3.4). The definition of the maximum deformation, $2R_{max}$, is somewhat subjective as it depends on identifying the first instance that the bag blows out behind the deformed droplet; therefore the empirical correlation of Zhao et al. (Reference Zhao, Liu, Li and Xu2010) (equation (3.4)) is used to be consistent with their work. All other inputs to (3.10) and (3.11) can be pragmatically identified and measured directly from experiments.

The characteristic breakup times for both the Rayleigh–Taylor and Rayleigh–Plateau instabilities are compared with the experiments in figures 4(a) and 4(b), respectively, where the experimental $t_b$ is measured as the time between the initiation of the rim ($T_i$, i.e. when the rim first forms and becomes susceptible to the instability) and the end of the breakup of the rim (i.e. when the instability fully penetrates and fragments the rim). Note that this definition is slightly complicated by the fact that it considers the breakup of the remaining rim between the nodes rather than of just the nodes themselves. This is necessary in the present analysis since the nodes do not typically fragment before the remaining rim. For this reason, the comparison of the instability breakup times is only approximate to test if the theoretical breakup times are within a reasonable range.

Figure 4. Comparison of the theoretical breakup time with the experimental breakup time for the (a) Rayleigh–Taylor ((3.9) given (3.10)) and (b) Rayleigh–Plateau ((3.9) given (3.11)) instability theories. The constant $C_t$ is fitted to the experiments for both cases.

The values of $C_t$ for both the Rayleigh–Taylor and Rayleigh–Plateau instabilities are obtained by fitting (3.9) (using (3.10) and (3.11) for the Rayleigh–Taylor and Rayleigh–Plateau cases, respectively) to the data, giving $C_t=7.1$ ($r^2=-0.04$) and $C_t=4.1$ ($r^2=0.55$) for the Rayleigh–Taylor and Rayleigh–Plateau instabilities, respectively. Both fitted values of $C_t$ are of the same order of magnitude as that given by Grant & Middleman (Reference Grant and Middleman1966), and also give $\zeta _0 = O($nm$)$, suggesting that both theories are reasonably plausible given the potential experimental error in $t_b$ and the increased perturbation due to the aerodynamic forces on the droplet and the rim. Although the Rayleigh–Taylor theory follows the data well for the BS morphology, the theory sharply over-predicts the breakup time through the bag morphology, resulting in $r^2<0$ for the whole fit. This likely arises due to the much smaller cross-sectional deformation that occurs in this regime (see (3.4)), which has a strong effect on the prediction of the acceleration, $a \propto (2R/d_0)^2$. However, we again stress that due to the challenges in estimating the experimental breakup time of the nodes compared with the rest of the rim, the comparison of the instability breakup times is only approximate and is meant to test if the theoretical breakup times are within a reasonable range.

3.2. Comparison of instability theories for rim node dynamics

Having shown that the Rayleigh–Plateau theory is equally as plausible as the well-used Rayleigh–Taylor theory for the instability of the rim that forms the nodes, we now turn to the prediction of the number of nodes formed (§ 3.2.1) and their size upon breakup (§ 3.2.2).

3.2.1. Node formation: number of nodes

Assuming that the toroidal rim is unstable to a wavelength $\lambda$, then the number of nodes, $N$, formed around the rim will be equal to the whole number of wavelengths that can fit around the circumference of the deformed droplet (Zhao et al. Reference Zhao, Liu, Li and Xu2010):

(3.12)\begin{equation} N = {\rm \pi}\frac{2R}{\lambda}. \end{equation}

For the case of the Rayleigh–Taylor instability, the substitution of (3.3) into (3.12) yields

(3.13)\begin{equation} N_{RT} = {\rm \pi}\frac{2R_{max}}{\lambda_{RT}} = \frac{\sqrt{C_D We}}{4} \left( \frac{2R_{max}}{d_0} \right)^2. \end{equation}

For the case of the Rayleigh–Plateau instability, the substitution of (3.5) into (3.12) yields

(3.14)\begin{equation} N_{RP} = {\rm \pi}\frac{2R_{max}}{\lambda_{RP}} = \frac{\rm \pi}{4.5} \frac{2R_{max}}{d_0} \frac{d_0}{h} \approx 0.69 AR, \end{equation}

where $AR=2R_{max}/h$ is the aspect ratio of the toroidal rim. Note that in this case, the maximum deformation, $2R_{max}$, has been used instead of the initiation deformation, $2R_i$. This is because while the instability is set up by the thickness of the rim at initiation, by the time the instability has grown enough to bifurcate the nodes of the rim, the rim has expanded closer to its maximum deformation. The rim thickness does not change significantly in the early stages of its expansion, and therefore the initiation thickness still governs the Rayleigh–Plateau instability (equation (3.5)) at the maximum deformation. Further detail is given in Appendix B.1.

Equations (3.13) and (3.14) are compared with the experimental measurements of the present analysis as well as the measurements of Dai & Faeth (Reference Dai and Faeth2001) and Zhao et al. (Reference Zhao, Liu, Li and Xu2010), which consist of water and ethanol droplet experiments, and the semi-empirical model of Zhao et al. (Reference Zhao, Liu, Li and Xu2010) in figure 5. Note that while the models are shown as continuous in $N$, they represent discrete steps in $N$, as $N$ must be a whole number due to the azimuthal symmetry condition of the toroidal rim.

Figure 5. Comparison of the present Rayleigh–Taylor (equation (3.13)) and Rayleigh–Plateau (equation (3.14)) models for number of nodes, $N$, versus $We$ compared with the present experiments for water droplets as well as the measurements of Dai & Faeth (Reference Dai and Faeth2001) and Zhao et al. (Reference Zhao, Liu, Li and Xu2010) for water and ethanol droplets and the semi-empirical model of Zhao et al. (Reference Zhao, Liu, Li and Xu2010). Filled and open markers indicate results from side- and end-on-view experiments, respectively.

Both models are found to follow the trend of the data and lie within its variability, with the Rayleigh–Taylor model (equation (3.13)) generally predicting the lower extent and the Rayleigh–Plateau model (equation (3.14)) generally predicting the mean. This suggests that the Rayleigh–Plateau instability may be more appropriate for modelling the instability of the rim than the Rayleigh–Taylor instability; however, considering the variation in $N$ for both the present data and that of Zhao et al. (Reference Zhao, Liu, Li and Xu2010) as well as the variation about the correlations for the extent of deformation ((3.4) and (3.8)), the results of figure 5 do not conclusively distinguish which mechanism is truly dominant. One of the possible causes for the experimental variation in $N$ is the coupling effect of the bag structures that form behind the droplet. For instance, in the BS morphology, the stamen can cause the bag to lose symmetry by inducing a fold in the bag (i.e. the beginnings of ‘twin-bag’ breakup; Jackiw & Ashgriz Reference Jackiw and Ashgriz2021), which may affect the symmetry of the rim, and thus its instability. Furthermore, the present Rayleigh–Taylor prediction, which uses the empirical correlation of Zhao et al. (Reference Zhao, Liu, Li and Xu2010) for $2R_{max}$, provides a better prediction of the trend than the Rayleigh–Taylor prediction of the same work, which instead uses the empirical correlation of Hsiang & Faeth (Reference Hsiang and Faeth1992). This highlights the importance of accurately modelling the extent of deformation on the downstream predictions such as $N$.

Note that the models presented here are simplified and neglect several potentially stabilizing effects. These effects are discussed further in Appendix B.1 and it is noted that these effects can only be compared with the present cases qualitatively.

3.2.2. Node breakup: node child droplet size

Once the nodes have formed on the rim of the droplet, they persist until the rim breaks, at which time they will break from the rim as child droplets having a size significantly larger than those formed by the rest of the rim. If all of the liquid contained within each wavelength on the liquid rim goes into the resulting droplets, then there would be no rim remaining after the node breakup. However, this is not what is observed, as a thinned region that undergoes its own breakup does remain between the nodes, likely due to nonlinear flow dynamics within the thin ligaments during the expansion of the rim. This is illustrated in figures 6(a) and 6(b). Therefore, it is only some fraction of the volume contained in each wavelength that contributes to the node child drops.

Figure 6. Illustration of node formation due to rim instability of wavelength $\lambda$ showing (a) the initial corrugation of the rim, (b) the formation of a node and (c) the assumed geometry of a cylindrical wavelength segment. The shaded areas denote the node volume.

By conservation of mass between a node droplet and the fraction $n$ of the cylindrical wavelength segment of the rim having diameter $h_i$ and length $\lambda$, as illustrated in figure 6(c), the rim child droplet size is found as

(3.15)\begin{equation} \frac{d_N}{d_0} = \left[\frac{3}{2} \left(\frac{h_i}{d_0}\right)^2 \frac{\lambda}{d_0} n \right]^{1/3}, \end{equation}

where $h_i$ is given by (3.6) and $\lambda$ is given by (3.3) and (3.5) for the Rayleigh–Taylor and Rayleigh–Plateau instabilities, respectively. For $n=1$, the entire wavelength segment will convert into the node child droplet. The problem, then, turns to determining how much of a given wavelength segment becomes part of the node (i.e. determining the appropriate value for $n$), and thus the node child droplet.

Since $n$ is the volume fraction of a wavelength segment of the rim that results in a node, $n$ must also be the fraction of the total volume of the rim that results in all of the nodes. Therefore, $V_N = n V_r$, where $V_N$ is the total volume of the node droplets and $V_r$ is the total rim volume. Assuming that the total volume of the rim does not change after its initiation, then one can use the expression provided by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021) for the volume contained in the rim at the initiation geometry as

(3.16)\begin{equation} \frac{V_{r}}{V_{0}}=\frac{3 {\rm \pi}}{2}\left[\left(\frac{d_i}{d_{0}}\right) \left(\frac{h_i}{d_{0}}\right)^{2}-\left(\frac{h_i}{d_{0}}\right)^{3}\right]. \end{equation}

Figure 7(a) plots the fraction of the total measured node volume (i.e. the sum of the node child droplet volumes) against the rim volume estimated by (3.16), where the dimensions of the rim are measured from experiments. The results show that under the assumption that the rim volume is constant after its initiation the node drops should account for roughly 90 % of the rim volume ($n=0.9$). However, this result is not physically accurate, as the remainder of the rim does not thin sufficiently to only account for 10 % of the total rim volume. Furthermore, the result $V_N/V_r>1$ is not physical if it is assumed that there is no flow between the inner portion of the disk and the rim. This suggests that liquid continues to feed the rim from the windward disk after the rim's initiation. Essentially, it is suggested that as the nodes draw liquid from the remaining rim, the remaining rim in turn draws liquid from the inner portion of the windward disk in order to maintain its equilibrium shape for as long as possible. As such, the rim scavenges liquid from the windward disk, and thus the total rim volume is better represented by the total volume of the windward disk, $V_d$. These flows are illustrated in figure 8. Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021) give an expression for the volume of the windward disk as

(3.17)\begin{equation} \frac{V_{d}}{V_0} = \frac{3}{2} \left[ \left(\frac{d_i}{d_0}\right)^2 \left(\frac{h_i}{d_0}\right) - 2 \left(1-\frac{\rm \pi}{4}\right) \left(\frac{d_i}{d_0}\right) \left(\frac{h_i}{d_0}\right)^2 \right]. \end{equation}

Figure 7. Node volume versus $We$ relative to (a) the predicted volume of the rim (equation (3.16)) and (b) the predicted volume of windward disk (equation (3.17)).

Figure 8. Illustration of the flows in the rim formation process, where liquid flows from the centre of the windward disk into the rim and from the rim into the disk.

Figure 7(b) shows that the mean volume of the nodes relative to the windward disk is $V_N/V_d\approx 0.4$, while the minimum volume is $V_N/V_d\approx 0.2$. As given earlier, the maximum physical value is $V_N/V_d\approx 1$. Assuming that a majority of the liquid in the windward disk will eventually contribute to the rim, then the volume fraction of the nodes relative to the windward disk will be representative of the volume fraction of the wavelength segment that becomes the node, $V_N/V_d \approx n$. Therefore, the minimum, mean and maximum values of $n$ are $n\approx 0.2$, $n\approx 0.4$ and $n=1$, respectively.

Using these values of $n$ and assuming the Rayleigh–Taylor (equation (3.3)) and Rayleigh–Plateau (equation (3.5)) wavelengths, (3.15) is compared with the experimental data in figures 9(a) and 9(b), respectively.

Figure 9. Node child droplet size versus $We$ compared to predictions of the (a) Rayleigh–Taylor (equation (3.3)) and (b) Rayleigh–Plateau (equation (3.5)) mechanisms.

Figure 9 shows that both the Rayleigh–Taylor and Rayleigh–Plateau instability theories with (3.15) give reasonable prediction of the range of node droplet sizes in the present experiments, where the mean node volume fraction, $n=0.4$, predicts the mean trend of the node breakup sizes, $d_N$, and the maximum and minimum values of $n=1$ and $n=0.2$, respectively, capture the upper and lower limits of the spread of $d_N$. Firstly, this indicates that the present analyses well represent the node formation and breakup process. Secondly, the variation in the node volume fraction, $n$, is a good indicator of the variation in sizes produced.

The predictions of the Rayleigh–Taylor and Rayleigh–Plateau instability theories are very similar. The Rayleigh–Taylor theory predicts $d_N$ somewhat better than the Rayleigh–Plateau theory, which slightly under-predicts the mean trend; however, the difference may be explained by the instability taking hold slightly after the rim is initiated, which would slightly reduce $h_i$ (3.6) owing to the thinning of the rim, and thus reduce the prediction of $d_N$. Identifying the instant that the instability takes hold, in particular in the time between the initiation ($i$) and the maximum (${max}$) extents of deformation, is not possible from the present experiments as the rim is not resolved with enough detail and dimension to quantify the development of its corrugation in the early stages of the instability. The results of figure 9 are therefore not sufficient to conclude which instability is truly dominant; however, at this time, the Rayleigh–Taylor instability theory appears to give the most accurate results for $d_N$ in the present experiments.

It is of note that the difference in the prediction of the Rayleigh–Taylor and Rayleigh–Plateau theories for $N$ is greater than for $d_N$ (compare figures 5 and 9). Although these two quantities could be well correlated if $d_N$ were derived from volume conservation for the whole rim, the derivation used here only considers volume conservation in one wavelength segment of the rim. As a result, $N$ and $d_N$ are not directly related to each other; however, they are related indirectly through the rim geometry parameters, $2R$ and $h$ (note that $\lambda$ is itself dependent on $2R$ and $h$). Comparing the predictions of $N$ ((3.13) and (3.14) for the Rayleigh–Taylor and Rayleigh–Plateau theories, respectively) with the predictions of $d_N$ (equation (3.15)), it can be seen that $d_N$ is less sensitive than $N$ to $2R$. Since the Rayleigh–Taylor and Rayleigh–Plateau theories are each better suited to a different definition of $2R$, the difference between the models for $N$ is greater than for $d_N$.

3.3. Remarks on the rim node instability

In the literature, the only instability that is considered in the growth of the rim nodes is the Rayleigh–Taylor instability (Zhao et al. Reference Zhao, Liu, Li and Xu2010; Xu, Wang & Che Reference Xu, Wang and Che2020; Qian et al. Reference Qian, Zhong, Zhu and Lin2021); however, in the present work, it was found that the Rayleigh–Plateau instability theory gives comparable results with those of the Rayleigh–Taylor instability in terms of breakup time, number of nodes formed and the resulting breakup sizes of the node drops. In terms of the breakup time and the number of nodes formed, the Rayleigh–Plateau theory was found to provide a slightly better prediction than the Rayleigh–Taylor theory, while for the final breakup size of the nodes, the Rayleigh–Taylor theory provided the better agreement with the experiments. In all cases, however, the variations were small enough that they were easily attributed to experimental errors or to slight variations in the definition of the extent of deformation at which the instability takes hold. As a result, we cannot conclude quantitatively which mechanism is truly occurring in the case of the node formation in the rim. Furthermore, other analyses of the instabilities of rims, such as the rim on the edge of a receding sheet, have also found that the two mechanisms provide similar results in predicting the instability of the rim (Lhuissier & Villermaux Reference Lhuissier and Villermaux2011; Wang et al. Reference Wang, Dandekar, Bustos, Poulain and Bourouiba2018). In fact, Lhuissier & Villermaux (Reference Lhuissier and Villermaux2011) concluded that the predictions from these two idealized mechanisms are ‘intrinsically undistinguishable’.

Although the results of these two theories may be functionally interchangeable for the present case in terms of a quantitative analysis, it may be instructive for future study to consider the qualitative differences between the two theories. In particular, a potentially severe flaw with the Rayleigh–Taylor instability approach is that it considers no changes to the original derivation for the case of the instability of an infinite planar interface when applied to the finite cylindrical interface of a rim. For the cylindrical rim, the additional curvature of the surface in the case of a cylindrical (or at least approximately cylindrical) body, in particular the axial curvature, would be expected to have a significant effect on the planar assumption of the applied Rayleigh–Taylor theory. Krechetnikov (Reference Krechetnikov2010) showed that the interplay of both the Rayleigh–Taylor and Rayleigh–Plateau mechanisms must be considered in the problem of the instability of a liquid sheet edge. Since the surface tension component stabilizes the Rayleigh–Taylor mechanism while destabilizing the Rayleigh–Plateau mechanism, as the scale of the problem decreases, the Rayleigh–Plateau mechanism becomes dominant. For the case of the rim at the edge of an expanding liquid sheet (i.e. in droplet impact), Wang et al. (Reference Wang, Dandekar, Bustos, Poulain and Bourouiba2018) found that, in most cases, the combined dynamics collapse to those of the Rayleigh–Plateau instability. Although the case of Wang et al. (Reference Wang, Dandekar, Bustos, Poulain and Bourouiba2018) is not directly transferable to the present one, it is illustrative of the effect of the curvature of the rim on the planar assumption of the Rayleigh–Taylor analysis.

In spite of the potential phenomenological downfall of the Rayleigh–Taylor theory described above, it is nonetheless unavoidable that, in the present work, the Rayleigh–Taylor theory gives a somewhat better prediction of $d_N$. Since the present work aims to develop a relationship for the distribution of droplet sizes originating from the variety of mechanisms of child droplet formation in aerodynamic droplet breakup, we will carry forward the Rayleigh–Taylor model; however, we stress that future works should not overlook the viability of the Rayleigh–Plateau theory for this case, as this mechanism may be more generalizable than the Rayleigh–Taylor theory if sufficient improvement in the modelling of the dynamics can be made.

3.4. Remaining rim dimensions

Since the breakup of the remaining rim will depend on its dimensions, it is necessary to determine these dimensions, in particular the thickness. Two factors contribute to the thinning of the rim. Firstly, since the rim major diameter expands in time, it must thin in order to conserve its volume. This is described in Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021), where the rim is found to thin as $h_f/h_i = \sqrt {R_i/R_f}$, where $R_f/R_i$ is the rim expansion ratio through the breakup phase. The second factor is the scavenging of the rim's volume by the formation of the nodes. Since the nodes form from the volume contained in the rim, the remaining rim segments would be expected to have a smaller volume, and thus a smaller thickness, than if the nodes had not formed. The initial and final volumes of the rim are related by the mass fraction remaining in the rim after the nodes have formed, $(1-n)$, as $V_{r,i} (1-n) \approx V_{r,f}$, where the rim volume, for either initial or final geometries, is $V_r = {\rm \pi}^2 (2R) h^2 / 4$.

However, this assumes that the volume of the rim, including the nodes, is constant after its initiation. In the previous section, it was shown that there is some evidence to suggest that liquid continues to flow into the rim after its initiation. We then define an additional parameter, $m$, that describes the volume fraction of the rim scavenged from the interior portion of the windward disk after the rim's initiation. The mass fraction remaining in the rim accounting for the flow from the interior portion of the windward disk after the nodes have formed is thus $(1-n+m)$, and $V_{r,i} (1-n+m) \approx V_{r,f}$. Rearranging and normalizing by the initial droplet volume gives

(3.18)\begin{equation} \frac{h_f}{d_0} = \frac{h_i}{d_0} \sqrt{(1-n+m) \frac{2R_i}{d_0} \frac{d_0}{2R_f}}. \end{equation}

For $n>m$, more liquid is flowing from the rim to the nodes than is flowing from the windward disk to the rim, while for $n< m$, more liquid is flowing from the windward disk to the rim than from the rim to the nodes. The case of $n=m$ is compared with the measurements in figure 10.

Figure 10. Measured final rim thickness, $h_f$, compared to the prediction of (3.18) from measurements of $h_i$ and $2R_i$ assuming $n=m=0.4$.

Figure 10 shows that the results of (3.18) for $n=m$ agree well with the present experimental measurements with no obvious dependency on $We$ for the studied range of $9< We<30$. An increase in $n$ relative to $m$ would result in a lower theoretical prediction, and thus worse agreement with the measurements. Therefore, the earlier suggestion that the rim draws additional liquid from the core of the windward disk after its initiation is supported. Conversely, a small increase in $m$ relative to $n$ may improve the present prediction, suggesting that there is a greater flow from the windward disk into the rim than from the rim into the nodes. Although this suggests an interesting finding in the flow dynamics of droplet breakup, the sensitivity of (3.18) to $n$ and $m$ is relatively small ($O(\sqrt {1-n+m})$), and thus, for simplicity in the remainder of the present analysis, the value of $n=m$ is used, recovering the result of Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021) despite the formation of the nodes, which were shown in § 3.2.2 to scavenge $n\approx 0.4$ of the rim's mass.

4.1. Liquid sheet recession and thickness

When a liquid sheet ruptures, the curvature of the sheet at the edge of the rupture results in a high surface tension force away from the rupture site, causing the edge of the sheet to recede away from the rupture site. As it does so, the edge of the sheet picks up more liquid from the sheet, forming a rim thicker than the sheet. Culick (Reference Culick1960) showed that these dynamics result in a constant retraction speed, $u_{rr}$, that depends on the sheet thickness, $h$, as

(4.1)\begin{equation} h = \frac{2 \sigma}{\rho_l u_{rr}^2}, \end{equation}

which is referred to as the Taylor–Culick law. Since the liquid sheet thickness is exceptionally difficult to measure, the Taylor–Culick law is often used to infer the sheet thickness based on the more easily measured retraction speed (Opfer et al. Reference Opfer, Roisman, Venzmer, Klostermann and Tropea2014; Poulain, Villermaux & Bourouiba Reference Poulain, Villermaux and Bourouiba2018; Villermaux Reference Villermaux2020).

In the present experiments, $u_{rr}$ is found by measuring the displacement of the edge of the receding rim relative to a relatively stationary point on the deformed droplet, such as the main rim, and taking its derivative. The parallax effect of the curvature of the bag tip is ignored for simplicity. To eliminate the effect of noise on the derivative, the displacement measurement is smoothed using a Savitzky–Golay filter of order 2 for a window length equal to the nearest odd number of measurement points, rounded down. By using this filter, it is assumed that the receding rim experiences a constant acceleration as it recedes, which matches the experiments well as shown in figure 12(a). An estimate of the error in the measurement is given by the error bars in figure 12(a), assuming that there is an error of $\pm 2$ pixels in identifying both the location of the sheet edge and the reference point (which itself is moving with the droplet), leading to a total error of $\pm 4$ pixels (about 0.17 mm). The bag recession speed is found to decrease linearly with time, indicating based on the Taylor–Culick law (equation (4.1)) that the bag is thinnest at its rupture site and thickest near the main rim of the deformed droplet, as shown in figure 12(b).

Figure 12. (a) Plot of the bag recession displacement in time, comparing experimental measurement with smoothed result using the Savitsky–Golay filter, and the recession speed calculated based on the derivative of the smoothed displacement. (b) Estimation of bag thickness based on calculated bag recession speed and Taylor–Culick law (equation (4.1)). Example case for a droplet having $d_0=2.07$ mm and $We=10.1$. Rupture ($t=0$ ms reference) at $T=2.08$ (6.95 ms).

The estimated bag thickness results are presented for the present experiments in figure 13, where the markers indicate the minimum thickness and the bars indicate the range of the thickness to its maximum.

Figure 13. Estimated sheet thickness versus $We$. Vertical bars show the time variation in $h$, where $h_{min}$ is indicated by the markers. The shaded region gives the standard deviation of $h_{min}$ across all experiments.

The present results indicate that there is no obvious relationship between $We$ and $h_{min}$, indicating that, at least for the present range of $We<30$ for water, there is no dependency of $h_{min}$ on $We$. Opfer et al. (Reference Opfer, Roisman, Venzmer, Klostermann and Tropea2014) found $h_{min}$ to be of a similar order of magnitude ($h_{min} = 1$–$0.1\ \mathrm {\mu }$m over the range $We=10$–$15$) despite the initial droplet sizes being approximately $>50\,\%$ larger ($d_0 = 3.2$ mm) than in the present study. It is therefore reasonable to conclude, given the present results, that there is a universal critical bag thickness at which the bag ruptures for water. Based on figure 13, the average of the minimum bag thickness is $\bar {h}_{min} = 2.3 \pm 1.2\ \mathrm {\mu }$m.

Opfer et al. (Reference Opfer, Roisman, Venzmer, Klostermann and Tropea2014) suggested an alternative approach to inferring the evolution of the bag thickness using a conservation of volume argument on the deforming bag where the bag size in the downstream direction at the time of burst is determined empirically. Their results indicated that the bag thickness decreases with increasing $We$ from $h_{min} = 1$ to $0.1\ \mathrm {\mu }$m over the range $We=10$ to $15$ and was found to be independent of the liquid properties over their tested range ($\mu _l = 1$–$16.1\ {\rm mPa}\,{\rm s}, \sigma =21.2\text {--}72$ mN m$^{-1}$). The $We$ dependency of the bag thickness is inherited from the $We$ dependency of the bag size at burst, which increases with $We$ over their tested range. While their prediction of the evolution of the bag thickness coincides with corresponding measurements using the Taylor–Culick sheet recession method, it assumes that mass is conserved in the bag throughout its deformation, to which evidence to the contrary was presented in section § 3.4. While the trend in the bag thickness results of Opfer et al. (Reference Opfer, Roisman, Venzmer, Klostermann and Tropea2014) is contrary to the constant value found in the present work, the order of magnitude and spread of their data is consistent with the present results. Since the actual mechanisms leading to the bag's rupture, and thus its terminal thickness, are still unknown, the present analysis will continue with the constant minimum bag thickness found earlier; however, the work of Opfer et al. (Reference Opfer, Roisman, Venzmer, Klostermann and Tropea2014) highlights the need for further study of the mechanisms of the bag's rupture.

The dynamics governing the rupture of surface bubbles considers, mainly, the competing effects of evaporation (which occurs at times scales $>10$ s) and film drainage at the meniscus of the bubble–surface interface (Bourouiba Reference Bourouiba2021). In some cases, Marangoni flows induced near the meniscus by variations in surfactant concentration or temperature can also affect the lifetime of the bubble (Poulain et al. Reference Poulain, Villermaux and Bourouiba2018). However, all of these effects have been studied at the relatively long time scales of surface bubble lifetimes ($O(1{\rm s})$) as compared with the lifetime of the bag in aerodynamic droplet breakup ($O({\rm ms})$) and do not account for the rapid stretching of the bag. Previous analyses of aerodynamic bag breakup have assumed that the bag ruptures when it develops a Rayleigh–Taylor instability at its tip that results in the piercing of the bag (Vledouts et al. Reference Vledouts, Quinard, Vandenberghe and Villermaux2016; Jackiw & Ashgriz Reference Jackiw and Ashgriz2021); however, this theory depends on $We$ and over-predicts the present results by an order of magnitude.

4.2. Instability of the receding rim

Having a view of the bag thickness and recession characteristics, the instability of the receding rim can now be studied. Here, we compare the universal rim thickness criterion of Wang et al. (Reference Wang, Dandekar, Bustos, Poulain and Bourouiba2018), previously shown to be applicable for both the edges of circular planar sheets and the rupture of surface bubbles, with the present case of aerodynamic droplet breakup. In this theory, the thickness of the receding rim, $b$, is given, universally, by the criterion $Bo=\rho _l b^2 a / \sigma = 1$. The acceleration used for the case of the surface bubble (i.e. a spherical cap) is taken as the centripetal acceleration:

(4.2)\begin{equation} a_c = \frac{u_{rr}^2}{R}, \end{equation}

where $R$ is the radius of the bubble. For the present case of the bag, $R$ is approximated by the bag size in the streamwise direction. The receding rim thickness is thus found, by rearranging $Bo=1$, as

(4.3)\begin{equation} b = \sqrt{\frac{\sigma}{\rho_l a_c}}. \end{equation}

The receding rim, having thickness $b$, is then unstable to the Rayleigh–Plateau instability, as

(4.4)\begin{equation} \lambda_{RP} = 4.5 b. \end{equation}

An important distinction in this model is that the centripetal acceleration sets the thickness of the receding rim but does not govern its instability. Several other theories have been proposed for the instability of the receding rim; however, we have found that the theory of Wang et al. (Reference Wang, Dandekar, Bustos, Poulain and Bourouiba2018) best matches the present results. An overview and comparison of several theories for the instability of receding rims are given in Appendix C.

When comparing instantaneous measurements of the receding rim wavelength, it is necessary to consider that the rim's geometry is changing in time; mainly, that the rim is expanding. As a result, the nodes may be measured to be farther apart than the estimated wavelength. However, since the rim is assumed to have a constant thickness, it will also continue to be unstable in the same sense throughout its expansion, allowing new nodes to form when the separation between two existing nodes approaches $2\lambda$. This establishes an upper limit on the measurement as $2\lambda$. The phenomenon of the nascent node is shown in figure 14. The nodes on the receding rim first form in figure 14(a) and begin to digitate in figure 14(b). In figure 14(c,d), the nodes move apart owing to the expansion of the receding rim. By figure 14(e), the nodes are far enough apart for a new node to form on the rim. In figure 14(f), the nascent node is clearly visible.

Figure 14. Nascent nodes on receding rim of a breaking water droplet. First, the primary nodes (a) form and (b) digitate. As the nodes move apart (c–e), they open up new space for a nascent node to form on the receding rim (f). The initial droplet has diameter $d_0=2.08$ mm at $We=12.2$. The first panel is at $T=2.18$ (6.95 ms) and the time separation between the panels is ${\rm \Delta} T=0.016$ (0.05 ms). See supplementary movie 2.

To predict the wavelength based only on the initial droplet conditions, the minimum bag thickness, $h_{min}$, and the bag radius, $R$, must be determined.

As mentioned in § 4.1, the mean measured minimum bag thickness is $\bar {h}_{min} = 2.3 \pm 1.2\ \mathrm {\mu }$m, and it is proposed that this is indicative of a universal minimum bag thickness for aerodynamic droplet breakup for water. Since the mechanism that leads to this minimum thickness is not known, the effect of viscosity cannot be determined without further investigation. Based on studies of the lifetime of surface bubbles (Poulain et al. Reference Poulain, Villermaux and Bourouiba2018), one can theorize that the minimum thickness could be either larger or smaller with higher viscosity. On the one hand, the stabilizing effect of higher viscosity has been shown to increase the lifetime of surface bubbles, which may permit more stretching and thus a thinner sheet. On the other hand, higher viscosity also impedes drainage in surface bubbles, leading to a thicker sheet (which is ultimately one of the stabilizing effects mentioned previously). It should also be noted that a higher viscosity will tend to stabilize the receding rim, and could possibly impede it from corrugating as described here.

The bag radius at burst, $R$, used to calculate the centripetal acceleration can be estimated based on the semi-analytical bag growth model of Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021). This model is based on a force balance between the thinning bag and the stagnation point flow at its centre, where the breakup point is predicted using the Rayleigh–Taylor piercing criterion of Vledouts et al. (Reference Vledouts, Quinard, Vandenberghe and Villermaux2016). From the analysis, Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021) find the bag size, $\beta$, as

(4.5)\begin{align} \frac{\beta(t^{*})}{d_{0}}&=\frac{3}{4}\left(\frac{V_{0}}{V_{b}}\right) \frac{1}{\tau^{2}}\left\{\left[\left(\frac{2 R_{i}}{d_{0}}\right)-2\left(\frac{h_{i}}{d_{0}}\right)\right]^{2} \frac{t^{* 2}}{2}\right.\nonumber\\ &\quad +\left.\left(\frac{2 \dot{R}}{d_{0}}\right)\left[\left(\frac{2 R_{i}}{d_{0}}\right)-2\left(\frac{h_{i}}{d_{0}}\right)\right] \frac{t^{* 3}}{3}+\left(\frac{2 \dot{R}}{d_{0}}\right)^{2} \frac{t^{* 4}}{12}\right\}, \end{align}

where $h_i$ and $d_i$ are given by (3.6) and (3.8), respectively, and the bag volume $V_b/V_0$ is estimated as

(4.6)\begin{equation} \frac{V_{b}}{V_{0}}=\frac{3}{2}\left[\left(\frac{2 R_i}{d_{0}}\right)^{2}\left(\frac{h_i}{d_{0}}\right)-\left(\frac{2 R_i}{d_{0}}\right)\left(\frac{h_i}{d_{0}}\right)^{2} \left(2+\frac{\rm \pi}{2}\right)+{\rm \pi}\left(\frac{h_i}{d_{0}}\right)^{3}\right]. \end{equation}

Equation (4.5) is evaluated at the burst time, $t^{\ast }_b$:

(4.7)\begin{equation} t^{{\ast}}_{b}=\frac{\left[\left(\dfrac{2 R_{i}}{d_{0}}\right)-2\left(\dfrac{h_{i}}{d_{0}}\right)\right]}{\left(\dfrac{2 \dot{R}}{d_{0}}\right)}\left[{-}1+\sqrt{1+C \frac{8 \tau}{\sqrt{3 We}} \sqrt{\frac{V_{b}}{V_{0}}} \frac{\left(\dfrac{2 \dot{R}}{d_{0}}\right)}{\left[\left(\dfrac{2 R_{i}}{d_{0}}\right)-2\left(\dfrac{h_{i}}{d_{0}}\right)\right]}}\right], \end{equation}

to find the bag size at burst used in the present analysis for the final radius of the bag, $R_f = \beta (t^{\ast })$. The value of $C$ is determined to be $9.4$, considering the modifications made in the present analysis to the prediction of $2R_i$ (see Appendix A). While this model neglects the flow of mass from the interior of the windward disk into the rim described in § 3.2.2, Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021) find that it is sufficient in predicting the size of the bag at burst. This model for $R_f$ is sufficient for $We<80$ where the bag grows in the downstream direction. However, for $We>80$ where the bag periphery is dragged downstream, a value of $R_f = R_i$ may be more suitable, as the bag rim does not significantly increase in diameter after initiation.

The model of Wang et al. (Reference Wang, Dandekar, Bustos, Poulain and Bourouiba2018) ((4.4) given (4.3)) using the empirical value of $\bar {h}_{min} = 2.3 \pm 1.2\ \mathrm {\mu }$m and the semi-analytical model $R_f = \beta (t^{\ast })$ ((4.5) given (4.7)) is compared with the present measurements in figure 15, where the solid line gives the model prediction given the mean value $\bar {h}_{min}=2.3\ \mathrm {\mu }$m and the grey area gives the limits owing to the standard deviation of $\pm 1.2\ \mathrm {\mu }$m. The dashed line gives the upper limit of $2\lambda$, owing to the nascent nodes on the expanding rim.

Figure 15. Comparison of the measured receding rim wavelength with the model of Wang et al. (Reference Wang, Dandekar, Bustos, Poulain and Bourouiba2018) ((4.4) given (4.3)) using the empirical value of $\bar {h}_{min} = 2.3 \pm 1.2\ \mathrm {\mu }$m and the semi-analytical model $\beta (t^{\ast }_b)=R$ (equation (4.5); Jackiw & Ashgriz Reference Jackiw and Ashgriz2021) with the present measurements. The solid line gives the model prediction given the mean value $\bar {h}_{min}=2.3\ \mathrm {\mu }$m, while the grey area gives the limits owing to the standard deviation of $\pm 1.2\ \mathrm {\mu }$m. The dashed line gives the upper limit of $2\lambda$, owing to the nascent nodes on the expanding rim.

The results of figure 15 indicate that the model of Wang et al. (Reference Wang, Dandekar, Bustos, Poulain and Bourouiba2018) ((4.4) given (4.3)), using the presently determined empirical value of $\bar {h}_{min} = 2.3 \pm 1.2\ \mathrm {\mu }$m and the semi-analytical model $\beta (t^{\ast }_b)=R$ (equation (4.5); Jackiw & Ashgriz Reference Jackiw and Ashgriz2021) gives a good prediction of the receding rim wavelength. Notably, the variation in the minimum bag thickness appears to account almost exactly for the variation in the measured wavelength on the receding rim. This suggests that variations in the receding rim wavelength are the direct consequence of variations in the minimum bag thickness. Additionally, all measured values are well below the upper limit of $2\lambda$, predicted by the birth of new nodes on the expanding rim, with the exception of low $We$. Since the bag is the largest and has the fewest ruptures at low $We$, the recession time at low $We$ is the longest. Thus, the receding rim wavelength is most likely to expand towards the limit of $2\lambda$ at low $We$, which matches the results.

4.3. Liquid sheet breakup

Although the fragments resulting from the breakup of the liquid sheet are too small to accurately size from the present images ($O(\mathrm {\mu }{\rm m})$), it is instructive at this time to discuss the mechanisms of their breakup in the context of the bag dynamics described above. In § 6, the relationships discussed here will be compared with data that include sizes from the breakup of the bags.

As with the rupture of surface bubbles, the fragments from the bag come from two phenomena: the fingering of the receding rim (Lhuissier & Villermaux Reference Lhuissier and Villermaux2012; Wang & Bourouiba Reference Wang and Bourouiba2021) and the collision of two receding rims (Néel, Lhuissier & Villermaux Reference Néel, Lhuissier and Villermaux2020), as shown in figure 16. For both scenarios, previous studies suggest that the droplets from the breakup of the sheet will lie in a range spanning between the sheet thickness, $h$, and the Rayleigh–Plateau breakup size of the receding rim, $1.89b$, with the added potential for satellite droplets forming from the breakup of both cases. Satellite droplets are discussed further in § 5. This range agrees with the order of magnitude of the droplets observed in the present experiments.

Figure 16. Images of the two dominant bag breakup mechanisms for a water droplet: (a) the digitation of receding rim and (b) the collision of receding rims.

5.1. Comparison of rim breakup mechanisms

As described in the previous section, in some cases the local angle of recession, $\theta _{rr}$, can be so great that the collision is insufficient to dominate the breakup of the rim. In these scenarios, both the collision and Raleigh–Plateau mechanisms can simultaneously dominate the breakup in different regions of the rim. This is exemplified in figure 21, where two sides of the rim appear to be affected by each mechanism independently. In figure 21(a), side 1 experiences a direct collision from the receding rim, leading to its breakup. In contrast, in figure 21(b), $\theta _{rr}$ is very large and the receding rim travels parallel to side 2 such that they do not collide. As a result, the disturbance imparted to side 2 is not as significant as that imparted to side 1. Note that due to irregularities in the rim thickness, the rim thickness of each side is different, as shown in table 2.

Figure 21. Comparison of breakup by (a) the collision mechanism and (b) the Rayleigh–Plateau instability for a water droplet. Note that the separated panels on the right are at a larger time interval than those on the left (${\rm \Delta} T=0.038$ (0.1 ms) versus ${\rm \Delta} T=0.19$ (0.5 ms)). The initial droplet conditions are $d_0 = 2.02$ mm and $We = 16.6$. The start times of each panel are (a) $T=2.05$ (5.35 ms) and (b) $T=2.39$ (6.25 ms). See supplementary movie 4.

Table 2. Comparison of measured breakup times of the two sides of the rim with the theoretical breakup time due to the Rayleigh–Taylor instability. The initial droplet conditions are $d_0 = 2.02$ mm and $We = 16.6$.

The first indication that these mechanisms affect the breakup independently is the difference in the breakup time of each side. Table 2 gives the measured breakup times of sides 1 and 2 as well as the predicted breakup time of the Rayleigh–Plateau theory (Ashgriz Reference Ashgriz2011),

(5.5)\begin{equation} t_{b, RP} = \frac{13}{0.97 \sqrt{\dfrac{\sigma}{\rho_l h_f^3}}}, \end{equation}

and the ratio between the measured and predicted breakup times. From table 2, the measured breakup time of side 2 is three times that of side 1, which can be seen in figure 21 where side 2 (figure 21b) remains stable for much longer than side 1 (figure 21a) after the sheet has fully recessed. Notably, the measured breakup time of side 1 is less than that of side 2 despite side 1 being thicker than side 2, which would otherwise tend to increase the breakup time as predicted by the Rayleigh–Plateau theory. This observation emphasizes the effect of the collision mechanism where the strong addition of a specific disturbance significantly decreases the breakup time of the rim. These findings are consistent with the present theory that the impact of the corrugated rim accelerates the breakup of the main rim.

The mean child droplet sizes, $\bar {d}_r$, from each side of the rim as well as the predictions of the present collision theory (equation (5.2)) and the Rayleigh–Plateau breakup theory (equation (5.1)) along with comparisons of each with $\bar {d}_r$ are given in table 3. The results of table 3 indicate that the breakup sizes of side 1 are most accurately predicted by the present collision theory, while those of side 2 are most accurately predicted by the Rayleigh–Plateau theory.

Table 3. Comparison of measured breakup sizes of the two sides of the rim with the theoretical breakup sizes predicted by the collision (equation (5.2)) and Rayleigh–Plateau instability (equation (5.1)) breakup mechanisms. The initial droplet conditions are $d_0 = 2.02$ mm and $We = 16.6$.

The results of both breakup time and size for the two sides of the rim shown in the present case provide strong evidence for the existence and validity of the proposed receding rim collision breakup mechanism. Furthermore, the results show how both mechanisms can occur in the same breakup event, which contributes to a wider distribution of sizes generated by the breakup.

Note that while the predicted breakup time of side 2 (the Rayleigh–Plateau case) is also larger than that measured, this is likely due to the additional disturbance imparted to the rim by the collision of the receding rim, which would result in a lower breakup time than predicted by (5.5); however, in this case, the additional disturbance only decreases the breakup time and does not significantly affect the wavelength that breaks the rim. This is evidenced by the breakup sizes from the rim.

5.2. Summary of rim breakup

Having a view of the various mechanisms of rim breakup, we can now give an overview of how the variation in mechanism results in the distribution of sizes that result from the breakup of the rim, which are shown in figure 22.

Figure 22. Range of rim child droplet sizes compared with the collision (with satellite droplets) and Rayleigh–Plateau mechanisms. The vertical lines denote the standard deviation of the measured rim droplet sizes for each case to show their spread. The grey area gives the range of the collision mechanism considering the standard deviation in the bag thickness (see § 4).

Firstly, the largest sizes in the rim's breakup result from the Rayleigh–Plateau mechanism. The prediction of (5.1) matches well the upper limit of the sizes in the present experiments in terms of both magnitude and trend, as shown in figure 22.

The majority of the breakup sizes result from the present collision mechanism, which well describes the mean size of the droplets. Two fundamental factors influence the collision mechanism sizes: the minimum bag thickness, $h_{min}$, which affects the receding rim instability, and the collision angle, $\theta _{rr}$, which affects the effective imposed wavelength. As discussed in § 4, $h_{min}$ is approximately constant for all present cases, as $h_{min} = 2.3 \pm 1.2\ \mathrm {\mu }$m, giving a small variation to the collision mechanism about the mean. This is denoted by the grey area in figure 22. The collision mechanism has the tendency to increase the breakup size as $\theta _{rr}$ increases.

The smallest sizes of the breakup come from the satellite droplets (equation (5.4)) that form from the head-on collision of the receding rim. Note that satellite droplets may also form from the oblique collision of the receding rim, as well as from the Rayleigh–Plateau rim breakup mechanism. These satellites will be larger than those from the head-on collision, and thus fill out the distribution of sizes between the smallest satellites and the primary breakup mechanisms.

6.1. Modelling of the breakup distributions

Following the previous works, it is assumed that the number of droplets formed by the mechanisms within each mode are related approximately by a two-parameter gamma probability density function, giving the number probability density function, $p_n$, as

(6.1)\begin{equation} p_n(x=d/d_0; \alpha, \beta) = \frac{x^{\alpha-1} \,{\rm e}^{{-}x/\beta}}{\beta^{\alpha} \varGamma(\alpha)}, \end{equation}

where $x=d/d_0$ are the droplet sizes, $\varGamma$ is the gamma function, $\alpha$ is the shape parameter and $\beta$ is the rate parameter. Both $\alpha$ and $\beta$ must be greater than zero. For the two-parameter distribution, $\alpha$ and $\beta$ can be estimated using the method of moments estimators as

(6.2)\begin{equation} \alpha = (\bar{x}/s )^2 \end{equation}

and

(6.3)\begin{equation} \beta = s^2 / \bar{x}, \end{equation}

where $\bar {x}$ and $s$ are the number mean and standard deviation of the distribution, respectively. The estimation of $\bar {x}$ and $s$ from the models presented in the present work is discussed later in this section.

Despite containing much less of the volume of the parent droplet, the smaller modes of breakup dominate $p_n$ and obscure the larger modes as the number of droplets resulting from the breakup mechanisms increases with decreasing child droplet size. It is therefore more helpful to compare the volume-weighted probability density functions, $p_v$, rather than $p_n$, as this allows significantly better visualization of the larger modes. Also, $p_v$ is more useful in most industrial processes as the primary functions of many processes involving sprays depend on how much of the sprayed volume breaks into a certain size range. Function $p_v$ is obtained by weighting $p_n$ by $x^3$ and renormalizing the distribution:

(6.4)\begin{equation} p_v = \frac{x^3 p_n(x)}{\displaystyle \int_{0}^{\infty} x^3 p_n(x) {{\rm d}\kern0.7pt x} } = \frac{x^3 p_n(x)}{\beta^3 \varGamma(\alpha+3) / \varGamma(\alpha)}. \end{equation}

The resulting distribution of the breakup will be the weighted sum of all modes of breakup, given by

(6.5)\begin{equation} p_{v,Total} = \sum_{i=1}^{\hat{n}} w_{i} p_{v,i}, \end{equation}

where $w_i$ is the volume weight of the $i$th mode of the breakup determined by conservation of volume and $\hat {n}$ is the total number of modes.

The modes that contribute to the breakup of the droplet depend on the breakup morphology. In particular, the existence and mechanisms of breakup of the undeformed core play a significant role in the overall breakup, as the core can contain a considerable volume of the parent droplet and undergoes a separate breakup event that contributes additional modes to the breakup (Jackiw & Ashgriz Reference Jackiw and Ashgriz2021). The methodology of modelling the multimodal distribution of aerodynamic droplet breakup (i.e. the estimation of $w_i$ for the modes of breakup as well as $\bar {x}$ and $s$ from the models presented in the present work) for the simpler case of the bag breakup morphology is outlined next, followed by the methodology for extending the model to the case of the breakup of the undeformed core in the MB breakup morphology.

6.2. Prediction of size distribution in bag breakup

In the case of the bag breakup morphology, no undeformed core is formed during the breakup. Consequently, the only modes that contribute to the breakup sizes in bag breakup are the rim nodes, the remaining rim and the bag sizes. The combined $p_v$ of the nodes, rim and bag is then given by the weighted summation of each mode:

(6.6)\begin{equation} p_{v,Total} = w_{N} p_{v,N} + w_{r} p_{v,r} + w_{b} p_{v,b}. \end{equation}

As discussed in § 3.2.2, the node volume is found as

(6.7)\begin{equation} w_N = \frac{V_N}{V_0} = \frac{V_N}{V_d} \frac{V_d}{V_0}, \end{equation}

where the mean value of $V_N/V_d \approx 0.4$ was determined in § 3.2.2 and $V_d/V_0$ is given by (3.17).

The volume of the rim,

(6.8)\begin{equation} w_{r} = V_r/V_0, \end{equation}

is given by (3.16). Note that, as discussed in § 3.4, the rim maintains an approximately constant volume despite the flow from the rim into the nodes owing to the flow into the rim from the bag. The volume of the bag is then given by the remaining volume as

(6.9)\begin{equation} \frac{V_b}{V_0} = 1 - \frac{V_N}{V_0} - \frac{V_r}{V_0}. \end{equation}

Although the the undeformed core is neglected in this instance, the definition of $w_N$, $w_r$ and $w_b$ still results in a leftover core volume of $V_c/V_0\approx 0.05$. Since this volume is neglected when determining the total distribution, the integral will be $\approx 0.95$. The justification for omitting this size from the present analysis is discussed later in this section.

As mentioned previously, the distribution parameters $\alpha$ and $\beta$ can be estimated from $\bar {x}$ and $s$ ((6.2) and (6.3), respectively), which, in turn, are related to the child droplet sizes of the breakup mechanisms presented in the present work.

In the previous sections, the sizes predicted for the mechanisms within each mode were shown to give reasonable estimations of the number-based mean and spread of the present data; however, these limits do not translate directly to the true mean ($\bar {x}$) and standard deviation ($s$) of the droplet sizes since the weighting of the mechanisms within each mode is not considered. This is evidenced by the distance from the estimated mean to the estimated upper and lower limits being asymmetric. The relative weighting is introduced by enforcing the two-parameter gamma distribution assumption (i.e. assuming that the mechanisms in each mode are weighted such that the mode overall follows a gamma distribution). This is done by assuming that the mean and standard deviation of the distribution are equal to the mean and standard deviation of all of the combined mechanisms (weighted equally) within each mode, which are modelled by discrete predictions. The mean and standard deviation are then given by

(6.10)\begin{equation} \bar{x} = \frac{\displaystyle \sum_{i=1}^{\check{n}} x_i}{\check{n}} \end{equation}

and

(6.11)\begin{equation} s = \sqrt{\frac{\displaystyle \sum_{i=1}^{\check{n}} (x_i - \bar{x})}{\check{n}}}, \end{equation}

where $x_i$ are the characteristic sizes of each mode ($x=d/d_0$) and $\check {n}$ is the total number of characteristic sizes.

The breakup mechanism of the node drops (§ 3) is the simplest of the three breakup modes of bag breakup as there is only one characteristic mechanism, which contains its own inherent variation. The variation in the Rayleigh–Taylor instability theory for the mechanism of the node breakup comes from the variation in the node-wavelength volume fraction, $n$, which ranges from its lower limit of $n=0.2$ to an upper limit of $n=1$ with a typical value of $n=0.4$. The number-based mean and standard deviation of the node mode are then estimated by taking the mean and standard deviation of the three characteristic sizes that result from each of these three values (i.e. the upper, lower and typical values) in (3.15), weighted equally. The three predicted characteristic breakup sizes of the nodes are compared with the distribution data of Guildenbecher et al. (Reference Guildenbecher, López-Rivera and Sojka2009) in figure 23(a).

Figure 23. Comparison of the present predictions for the multimodal distribution with measurements from Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017) for ethanol drops of initial diameter $d_0=2.54$ mm at $We=13.8$ for the breakup of (a) the nodes, (b) the remaining rim and (c) the bag.

The breakup of the rim (§ 5) is somewhat more complex, as it consists of two mechanisms, i.e Rayleigh–Plateau (equation (5.1)) and the receding rim collision (equation (5.2)), as well as the formation of satellite droplets and the inherent variation of each. However, the sensitivity of the mechanisms to their variations is relatively small, and thus it will be assumed that only the primary mechanisms, as well as the satellite droplets, dominate the number distribution. Therefore, the number-based mean and standard deviation for the breakup of the remaining rim are estimated from four sizes, weighted equally: the Rayleigh–Plateau mechanism ($d_{r,RP}$, (5.1)), the receding rim collision mechanism ($d_{r, coll.}$, (5.2), assuming $\theta _{rr}=0$) and the satellite sizes that result from each ($d_{r, RP\ sat.}$ and $d_{r, coll.\ sat.}$, (5.4)). The four predicted characteristic breakup sizes of the rim are compared with the distribution data of Guildenbecher et al. (Reference Guildenbecher, López-Rivera and Sojka2009) in figure 23(b).

The mechanisms of the bag's breakup (§ 4) are somewhat ill-defined, as the present experiments were unable to capture these sizes. The main mechanism identified in the literature is the Rayleigh–Plateau breakup of the receding rim ($d_{b,RP}=1.89b$, where $b$ is given by (4.3)), which will also be accompanied by satellite droplets ($d_{b,\ RP sat.}$, given by (5.4)). However, comparing with the data of Guildenbecher et al. (Reference Guildenbecher, López-Rivera and Sojka2009), the sizes may be much smaller and approach the thickness of the sheet, $h_{min}$. As such, we suggest that two additional mechanisms of breakup in the bag result in sizes of the order of the sheet thickness, $h_{min}$, as well as the rim thickness, $b$. The number-based mean and standard deviation of the bag mode are then estimated by taking the mean and standard deviation of these four characteristic sizes, weighted equally. The four predicted characteristic breakup sizes of the bag are compared with the distribution data of Guildenbecher et al. (Reference Guildenbecher, López-Rivera and Sojka2009) in figure 23(c).

Figure 23(a) shows that the largest peak is due to the node drops, where the characteristic sizes are given by $n=1$, $n=0.4$ and $n=0.2$. Figure 23(b) shows that the central peak is the result of the breakup of the remaining rim, where the characteristic sizes are given by the collision mechanism, the Rayleigh–Plateau instability and the satellite droplets resulting from both. Figure 23(c) shows that the smallest peak is due to the breakup of the bag where the characteristic sizes are given by the Rayleigh–Plateau breakup of the receding rim and its associated satellite droplets as well as the characteristic sizes of the receding rim and bag sheet thickness. For all three modes, the spread of each mode is well captured by the constituent breakup mechanisms.

Having estimated $\bar {x}$ and $s$ from the characteristic breakup mechanisms, $\alpha$ and $\beta$ can be determined using (6.2) and (6.3), respectively, and thus the combined $p_v$ of the three modes (equation (6.6)) can be predicted. Figure 24(b) compares the analytical prediction of the combined distribution with the distribution data of Guildenbecher et al. (Reference Guildenbecher, López-Rivera and Sojka2009). The distribution parameters determined by the present method are provided in table 4.

Figure 24. Comparison of the present predictions for the multimodal distribution with measurements from Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017) for ethanol drops of initial diameter $d_0=2.54$ mm at $We=13.8$. (a) Distribution predictions using analytical weighting. (b) Distribution with best-fit weighting.

Table 4. Distribution parameters for the case of ethanol drops of initial diameter $d_0=2.54$ mm at $We=13.8$, corresponding to the experiments of Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017) shown in figure 24.

From figure 24(a), it can be seen that while the location and width of the peaks are reasonably matched by the semi-analytical distribution, the heights of the peaks are not so well predicted. The contribution of the central peak due to the rim breakup is over-predicted, while the left and right peaks, due to the breakup of the bag and nodes, respectively, are under-predicted. One possible explanation for this is that the assumed distribution weights ((6.7)–(6.9)) do not accurately represent the experiment. The weights used in the present analysis were derived specifically for the bag breakup morphology, while the experiment is near the transition to the BS breakup morphology ($We \approx 15$; Jackiw & Ashgriz Reference Jackiw and Ashgriz2021). Near this transition, the breakup may exhibit the ‘twin-bag’ transitional morphology described in Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021), where a fold appears across the bag dividing it in two owing to the slightly larger mass at the centre of the bag due to the small undeformed core, which was neglected in the present analysis as discussed earlier. While this fold primarily affects the bag, it may also have an influence on the rim; in particular, the formation of additional nodes on the rim where it is met by the fold of the bags. The presence of this fold then changes the assumed dynamics of the rim and thus may reduce the remaining rim's volume and thickness more than presently predicted, leading to the over-prediction in the volume contribution of the rim. The fold itself may also contribute additional large sizes to the breakup as it may form as a ligament that is left behind after the rupture of the bag. The ligament left behind is visible in the digital in-line holography images of Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017), although it is not included in their breakup data, which are used here. This ligament may also be akin to the undeformed core neglected from the present analysis, which contains approximately 5 % of the initial droplet's volume corresponding to a size of $d/d_0\approx 0.38$. Furthermore, it may be that several of the experimental cases led to the formation of a small stamen that additionally contributes to the volume in the large mode. Another explanation for why the remaining rim volume weight is over-predicted is that the dynamics of the node formation varying at low $We$ is not fully captured. If the flow from the bag to the rim is insufficient to match the flow from the rim to the nodes (i.e. $m< n$ in (3.18); see § 3.4 and figure 8) then the volume contribution of the remaining rim will be lower while that of the nodes will be higher. This is the case at very low $We$ ($We<11$) where the present model gives a non-physical prediction for the bag volume of $w_b<0$ due to the over-prediction of the remaining rim volume, $w_r$. This suggests that there is a limiting factor in the flow from the bag to the rim at very low $We$ that is not captured in the present work; however, such an investigation is outside of the present scope.

To see the effect of having a better prediction of the distribution weights, $w_N$, $w_r$ and $w_b$ are fitted to the experimental data in figure 24(b). The fitted weights are given in table 4. Although the combined distribution naturally matches the experimental data better, there are some deviations such as the valleys between the peaks and the right tail of the large mode, which reaches above a value of $d_c/d_0=0.78$, giving drops that contain more than half of the initial droplet volume. Notably, in figure 9, although the prediction of the upper limit continues to increase as $We$ decreases, the upper bar of the lowest $We$ data points does not exceed $d_N/d_0\approx 0.7$. These discrepancies may be solved by a narrower distribution, which can be obtained by changing either of the gamma distribution parameters, $\alpha$ or $\beta$. Due to the nature of the probability density function, the height of the peaks may also be affected by these changes. However, it is important to remember that these distribution parameters are an inherited part of the core assumption of this section: that the relative contribution of each mechanism to the overall mode is approximated by a gamma distribution function. The proposition of altering either of these parameters thus must also include the proposition of selecting a different distribution entirely. In particular, a truncated or bounded distribution may be more appropriate, as some of the mechanisms presented in the present work offer an upper physical bound on the child droplet size that is not captured by a gamma distribution. However, the comparison of different distributions is left out of the scope of the present work.

6.3. Breakup of the undeformed core

At higher-$We$ droplet breakup, the undeformed core also contributes to the breakup of the parent droplet and may undergo its own multimodal breakup, i.e. the ‘core breakup’. Thus, the undeformed core contributes additional modes to the breakup, and must be included in the analysis. The volume of the undeformed core is given by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021) as

(6.12)\begin{equation} w_c = \frac{V_c}{V_0} = 1 - \frac{V_d}{V_0}, \end{equation}

where $V_d/V_0$ is given by (3.17). The bag volume is then the volume of the windward disk not contained in the nodes or the remaining rim:

(6.13)\begin{equation} w_{b} = \frac{V_b}{V_0} = \frac{V_d}{V_0} - \frac{V_N}{V_0} - \frac{V_r}{V_0}. \end{equation}

The volume contributions of the node and rim modes are calculated as in § 6.2 ((6.7) and (6.8), respectively). The first breakup event, wherein the initial windward disk breaks and leaves the undeformed core behind, is modelled as in § 6.2, while the core breakup is treated separately.

The predicted distribution of the first breakup consisting of the node, rim and bag breakup sizes is compared with the data of Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017) in figure 25(a). Note that the plotted distributions are weighted by their respective volume contributions with respect to the parent droplet; thus, their respective integrals are $<1$. The distribution parameters and weights are given in table 5.

Figure 25. Comparison of present predictions of the (a) first and (b) core breakup event measurements from Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017) for ethanol drops of initial diameter $d_0=2.55$ mm at $We=55.3$. Note that the distributions in (a) and (b) are weighted by their volume contributions relative to the parent droplet, and thus their respective integrals are $<1$.

Table 5. Distribution parameters for the case of ethanol drops of initial diameter $d_0=2.55$ mm at $We=55.3$, corresponding to the experiments of Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017) shown in figures 25 and 26.

Figure 25(a) shows that the first breakup contributes predominantly to smaller sizes. However, the contributing volume of the bag appears to be over-predicted, especially considering that the first breakup contributes less than half of the total breakup volume in the present case. This is likely due to the folds between the azimuthal bags being neglected. The folds between the azimuthal bags leave behind ligaments after the rupture of the bags, which result in sizes similar to the breakup of the rim. As a result, some of the mass of the bag breaks into sizes similar to those of the rim breakup. Since this effect is neglected, the smallest sizes due to the breakup of the bag are over-predicted, while the sizes near the characteristic sizes of the rim breakup are under-predicted.

After the first breakup event, the undeformed core remains, which may also break due to the aerodynamic forces. It is therefore necessary to analyse the breakup of the undeformed core in the same manner as the breakup of the parent droplet; however, assumptions must be made about the aerodynamic conditions of the undeformed core.

The undeformed core is assumed to be represented by an effective droplet size by

(6.14)\begin{equation} \frac{d_{c}}{d_0} = \left( \frac{V_c}{V_0} \right)^{1/3}, \end{equation}

where $V_c/V_0$ is given by (6.12). Since the core droplet does not form until after the first breakup event has occurred, the acceleration of the parent droplet must be considered. The speed of the core droplet after the first breakup event, $U_c$, is found by assuming a constant acceleration due to drag as

(6.15)\begin{equation} U_c = \frac{3}{4} \bar{C}_D T_b \left( \frac{2R_i}{d_0} \right)^2 \frac{d_0}{\tau}, \end{equation}

where $T_b$ is the non-dimensional burst time ($T_b = T_i + T^{\ast }_b$, where $T_i \approx 0.9$ (Jackiw & Ashgriz Reference Jackiw and Ashgriz2021) and $T^{\ast }_b$ is computed from (4.7)) and $\bar {C}_D$ is the average drag coefficient over the initiation and bag blow-out given by Chou et al. (Reference Chou and Faeth1997) as 3.33 (calculated from the weighted average of the reported drag coefficients for the initiation and bag blow-out times). Using $d_c$ and the relative velocity $U - U_c$, the conditions of the core droplet, $We$ and $\tau$, can be determined. The breakup model is then repeated until the core droplet does not meet the conditions for breakup ($We\approx 8.8$; Jackiw & Ashgriz Reference Jackiw and Ashgriz2021). Since the second breakup event constitutes only a portion of the initial parent droplet, the weights of the modes from the core breakup must be adjusted so that they represent the weighting with respect to the parent droplet by multiplying each of the core breakup weights by the core volume weight, $w_c$.

For the present case, the core droplet properties are found as $d_c = 2.15$ mm, $U_c = 7.8$ m s$^{-1}$ ($U-U_c = 13.2$ m s$^{-1}$) and $We = 18$. The present prediction of $U_c$ is in reasonable agreement with the speed of the larger child drops measured by Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017) ($U_c \approx 10$ m s$^{-1}$). The core droplet, therefore, undergoes BS breakup. In this case, $We$ of the core droplet is low enough that a second core breakup event does not occur; however, it is possible that a second core droplet that does not break will remain. The volume of the leftover undeformed core is implicit in the calculation of the weights of the core modes, resulting in a volume deficiency of 0.09 between the total core volume and the volume of the modelled modes. This is the volume of the leftover core droplet, which would have an estimated size of $d_c/d_0 \approx 0.45$. This second core droplet would have $We \approx 2$, and therefore would not break further. Such a droplet is visible in the images of Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017), but is not captured in the measurement with which the present analysis is being compared. For this reason, the second undeformed core is neglected in the present analysis.

The resulting size distribution of the core breakup event is compared with the results of Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017) in figure 25(b). As with figure 25(b), the plotted distributions are weighted by their respective volume contributions with respect to the parent droplet; thus, their respective integrals are $<1$. The distribution parameters and weights are given in table 5. Figure 25(b) shows how the core breakup contributes to the larger sizes in the size distribution of the breakup. The combined distribution of the first and core breakup events is compared with the data in figure 26(a).

Figure 26. Comparison of the present predictions for (a) the combined distribution and (b) the combined lumped distribution with measurements from Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017) for ethanol drops of initial diameter $d_0=2.55$ mm at $We=55.3$.

In the combined distribution, the largest sizes come from the nodes of the breakup of the core, while the smallest sizes come from the breakup of the bag in the first breakup event. The remaining modes fill in the intermediate sizes of the distribution, with reasonable agreement between the prediction and the data of Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017) for sizes above $d/d_0\approx 0.2$. However, the combined distribution suffers the same overestimation of the smallest sizes as the first breakup event as mentioned earlier in this section. With the correction of these factors, the predicted distribution would be expected to match the experimental data much better in this range.

Since the modes of the first nodes and second rim overlap, as well as the modes of the first rim and the first and second bags, it may be convenient to lump these modes together such that the overall breakup consists of only three modes instead of six: the first rim and the first and second bags giving the smallest mode, the first nodes and the second rim giving the intermediate mode and the second nodes giving the largest mode. In doing so, the effect of the folds between the azimuthal bags that results in the misprediction of the bag and rim contributions can be reduced. The distribution parameters of the combined bag and rim mode are then calculated by the weighted mean and standard deviation of the constituent breakup mechanisms of both the bag and rim modes of breakup. The weighted mean, $\bar {x}_w$, and standard deviation, $s_w$, are calculated as

(6.16)\begin{equation} \bar{x}_w=\frac{\displaystyle \sum_{i=1}^{\check{n}} w_{i} x_{i}}{\displaystyle \sum_{i=1}^{\check{n}} w_{i}} \end{equation}

and

(6.17)\begin{equation} s_w = \sqrt{\frac{\displaystyle \sum_{i=1}^{\check{n}} w_{i}\left(x_{i}-\bar{x}_{w}\right)^{2}}{\displaystyle \frac{(\check{n}-1)}{\check{n}} \sum_{i=1}^{\check{n}} w_{i}}}, \end{equation}

respectively, where the characteristic sizes of each mechanism are calculated as before, and the weights are calculated by dividing the mode weight by the number of mechanisms that contribute to that mode (i.e. each mechanism within a mode is weighted equally, as before). The combined lumped distributions are compared with the data of Guildenbecher et al. (Reference Guildenbecher, Gao, Chen and Sojka2017) in figure 26(b). As expected, the distribution that results from lumping overlapping modes together gives a better prediction of the distribution for $d/d_0 < 0.2$, as the effect of the misprediction of the contributions of the bag and rim modes is reduced.

6.4. Summary of distribution prediction

In the present work, the distribution of sizes resulting from droplet breakup is assumed to come from a combination of a plethora of mechanisms as well as subtle variations within each mechanism. These mechanisms are grouped in terms of the primary modes of the breakup (i.e. the geometries within the deformed droplet) such that the overall distribution is the result of the sum of the modes, which can lead to a multimodal distribution. While the volume weighting of each mode is determined by models for the volume estimation of each geometry, the relative number weighting of the mechanisms within each mode is assumed to vary as a two-parameter gamma distribution function where the parameters are estimated based on the mean and standard deviation of the predicted sizes for each mechanism within the mode. In bag breakup, there are three modes resulting from the rim nodes, the remaining rim and the bag. In MB breakup, there is the addition of the breakup of the undeformed core, which itself breaks into an additional three modes. However, in MB breakup, several of the modes overlap. When two or more modes overlap such that they appear as one, the weighted mean and standard deviation of the mechanisms contributing to each mode can be used to determine a lumped distribution that captures the overlapping modes as one. Figure 27 gives a flowchart of the calculation of the multimodal distribution.

Figure 27. Flowchart for the calculation of the multimodal size distributions.