3.1. Effect of nozzle geometry and scale

The design of the nozzle is one of the most important considerations in spray atomization as it governs the interaction between the liquid and air flows, which ultimately leads to the breakup of the liquid. Twin-fluid nozzle design can be difficult to analyse as one must consider both the gas and liquid flow geometries, as well as how they interact. As illustrated in figure 2, for concentric, annular twin-fluid nozzles, there are three primary dimensions: the liquid orifice diameter, $d_l$, and the inside and outside gas orifice diameters, $d_{g,i}$ and $d_{g,o}$, respectively. These dimensions dictate the liquid and gas velocities for a given mass flow rate or upstream pressure condition. Conventionally, it is assumed that the liquid jet will take on the diameter of the liquid orifice, $d_l$; however, Machicoane et al. (Reference Machicoane, Bothell, Li, Morgan, Heindel, Kastengren and Aliseda2019) and Kumar & Sahu (Reference Kumar and Sahu2020) found that at practical operating conditions the air flow over the step in the nozzle face between $d_{g,i}$ and $d_l$ results in a low pressure region that, with the aid of surface tension, causes the liquid jet to wet across the face of the nozzle such that its diameter will always be $d_{g,i}$. This is shown in figure 10, where images of the liquid jet without (main image) and with (overlaid ‘ghost’ image) gas flow are compared, showing how the liquid jet assumes the diameter $d_{g,i}$ when the gas flow is present.

Figure 10. Images showing liquid wetting across the nozzle face. The main image shows the case with no gas flow, where the liquid jet emerges at the liquid orifice diameter. The overlaid ‘ghost’ image shows the case with gas flow, where the liquid jet wets across the face of the liquid nozzle.

When the wetting effect occurs, the wetted diameter should therefore be considered over the liquid orifice diameter. This conclusion is supported by figure 11, where the spray measurements of the 2050-70 and 2850-70 nozzles, which are identical apart from $d_l$ (see table 1), are compared in terms of both (a) the spray $d_{32}$ for varying gas flow rate and (b) the size distribution for fixed gas flow rate at $\dot {m}_l=0.67$ g s$^{-1}$ and $x=80$ mm, showing functionally identical results. Therefore, when designing both nozzles and models for their atomization characteristics, $d_{g,i}$ should be taken as the important dimension governing the liquid jet rather than $d_l$. Note that this selection will also affect the estimation of the average liquid jet velocity as the flow area is larger, resulting in a lower average velocity than would have been estimated using $d_l$. However, it should also be mentioned that such an average velocity is not necessarily an accurate representation of the liquid interface velocity, as the jet will exhibit a more complex profile, especially in the wetting case. Nevertheless, the average liquid jet velocity serves as a useful characteristic measure of the liquid jet behaviour.

Figure 11. (a) Comparison of $d_{32}$ of sprays from 2850-70 and 2050-70 nozzles at varying $G_g$ and of (b) volume-frequency size histogram at fixed $G_g=319$ kg m$^{-2}$ s$^{-1}$. $G_l=526$ kg m$^{-2}$ s$^{-1}$ and $x=80$ mm for both plots. The arrow call-out in (a) indicates the conditions for which the size distributions are plotted in (b). The histograms are plotted in (b) as marked lines for clarity, where the height and centres of the bins are given by the markers.

Three cases arise as exceptions to this phenomenon. Firstly, in the extreme case of very low liquid flow and high gas flow, the liquid jet may be atomized faster than it is replenished, preventing the liquid flow from forming a true jet and resulting in other dynamics such as an intermittent spray. Secondly, if the difference between $d_l$ and $d_{g,i}$ is sufficiently large, then there may be no significant liquid wetting across the face of the nozzle and the liquid jet will maintain the diameter of the liquid orifice, $d_l$. Such is likely the case in Varga et al. (Reference Varga, Lasheras and Hopfinger2003) where no liquid wetting is obvious for the small liquid orifice nozzle ($d_{g,i}= 1.3$ mm with $d_l = 1$ mm and $0.32$ mm). Finally, in cases where the step thickness is much smaller than the orifice diameter, then the difference between $d_l$ and $d_{g,i}$ will be negligible, as is the case in Marmottant & Villermaux (Reference Marmottant and Villermaux2004b).

Since the liquid jet diameter is dependent on the gas inner jet diameter, the only way to practically alter the liquid jet diameter is to increase the overall size of the nozzle; however, changing the overall nozzle size will also affect its airflow, thus, the effect of the liquid and air flows must be considered in conjunction. The measured droplet size distributions for both the 2850-70 and 60100-120 nozzles for two different values of $G_g$ are compared for similar $G_g$ and $G_l$ in figure 12. Images of the spray near the nozzle for these cases are shown in figure 13. In both cases, the larger nozzle (60100-120, (c,d)) generates larger droplets than the smaller nozzle (2850-70, (a,b)). Additionally, the larger nozzle exhibits a bimodal droplet size distribution at lower $G_g$ than the smaller nozzle (compare figures 12a and 12c). Both of these effects can be attributed to the larger diameter liquid jet being more difficult to atomize than the smaller diameter liquid jet. When considering only the $d_{32}$ of the spray, the effect of the liquid diameter on the breakup would appear to be much more significant at lower $G_g$, as shown in figure 14. However, this is mainly due to the bimodal distribution that the larger nozzle produces over this range. In cases where the distribution is monomodal for both nozzles, the $d_{32}$ will be only slightly larger for the larger nozzle than for the smaller nozzle at the same $G_l$ and $G_g$. Additionally, it is important to note that for the extreme cases of high $G_l$ and low $G_g$, where the atomization is extremely poor and the larger mode dominates significantly over the smaller mode, the distributions measured by the Malvern Spraytec are clipped at 1000 $\mathrm {\mu }$m (e.g. figure 9e). As a result, the $d_{32}$ measurements provided by the instrument at such conditions will be underpredicted. While such conditions are included on the $d_{32}$ plots here for completeness, it is important to keep this limitation in mind for these extreme cases. The bimodal cases are highlighted in the plots of $d_{32}$ by open markers, where the affected cases are typically only the lowest two $G_g$ values at high $G_l$.

Figure 12. Volume-frequency droplet size histograms of the sprays from the 2850-70 (a,b) and 60100-120 (c,d) nozzles at two similar $G_g$ (a,c and b,d) for similar $G_l$. Conditions are indicated on the figure. For all cases, $x=80$ mm. Histograms correspond to images in figure 13.

Figure 13. Images of sprays from the (a,b) 2850-70 and (c,d) 60100-120 nozzles at two similar $G_g$ (a,c and b,d) for similar $G_l$. Conditions are indicated on the figure. For all cases, $x=80$ mm. Images correspond to the histograms in figure 12. Results are shown for (a) 2850-70: $G_g=319$, $G_l=526\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (b) 2850-70: $G_g=535$, $G_l=526\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (c) 60100-120: $G_g=328$, $G_l=493\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; and (d) 60100-120: $G_g=527$, $G_l=493\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$.

Figure 14. Plot of $d_{3,2}$ vs $G_g$ comparing the sprays of the 2850-70 and 60100-120 nozzles at $x=80$ mm for similar $G_l$ (indicated on the figure). Closed and open markers indicate the conditions at which the atomization is complete or incomplete, respectively.

While the presently reported trends in spray size with nozzle scale are consistent with established empirical correlations such as Nukiyama & Tanasawa (Reference Nukiyama and Tanasawa1950) and Rizk & Lefebvre (Reference Rizk and Lefebvre1983), Varga et al. (Reference Varga, Lasheras and Hopfinger2003) found that a smaller liquid orifice produces a mean droplet size slightly larger (O($\mathrm {\mu }{\rm m}$)) than the larger liquid orifice, which was attributed to a longer gas boundary layer attachment length in the nozzle with smaller $d_l$. In their work, Varga et al. (Reference Varga, Lasheras and Hopfinger2003) compared the $d_{32}$ of two nozzles with differently sized liquid jets at the same $u_g$ for the same $m_r$, i.e. the same liquid mass flow rate (note that while they named their mass-ratio term as the ‘mass flux ratio’, the actual ratio presented is the mass flow rate ratio, which does not normalize to the flow areas). The same comparison is shown for the present results in figure 15, which yields the same result: that the larger nozzle produces slightly smaller droplets on the order of a few $\mathrm {\mu }$m. However, rather than the boundary layer attachment length argument of Varga et al. (Reference Varga, Lasheras and Hopfinger2003), this result is attributable to the differences in $u_l$ of the two cases. By conservation of mass, a larger liquid jet operating at the same $\dot {m}_l$ will have a lower $u_l$ than a smaller jet, which, as discussed previously, tends to decrease the droplet sizes in the spray and opposes the effect of increasing the jet diameter.

Figure 15. Plot of $d_{32}$ vs $G_g$ comparing the sprays of the 2850-70 and 60100-120 nozzles at $x=80$ mm at $\dot {m}_l=1.67$ g s$^{-1}$. Closed and open markers indicate the conditions at which the atomization is complete or incomplete, respectively.

4.1. Initial interaction of the liquid jet with the coaxial gas stream

The first aspect that must be considered when modelling the atomization of a coaxial, twin-fluid spray is how the liquid jet initially interacts with the co-flowing air jet and how this interaction leads to the breakup of the liquid jet under the aerodynamic forcing of the gas flow. The two parameters that must be determined for modelling the breakup of the waves are the diameter of the liquid jet when it is exposed to the air stream, $d_{l}'$, and its relative velocity to the air stream, $u_r$.

Previous models assumed that the initial instability of the liquid jet was a leading term in this process. Due to the large difference in speeds between the liquid and gas flows, the liquid jet becomes susceptible to a Kelvin–Helmholtz shear instability. Previous works have shown that the finite-shear layer derivation of the Kelvin–Helmholtz instability for a planar interaction between a high-speed gas flow and a low-speed liquid flow is readily applicable to the configuration of a coaxial twin-fluid nozzle (Raynal Reference Raynal1997). From this basis, the Kelvin–Helmholtz wavelength, $\lambda _{KH}$, is derived as

(4.1)\begin{equation} \lambda_{KH} = 2 C_{\lambda} \frac{b_g}{\sqrt{Re_{bg}}} \sqrt{\frac{\rho_l}{\rho_g}}, \end{equation}

where $Re_{bg} = \rho _g u_r b_g / \mu _g$ is the Reynolds number of the annular gas flow based on the annular air-gap thickness, $b_g = d_{g,o} - d_{g,i}$, and $C_{\lambda }$ is a constant relating to the nozzle's contraction ratio, which has a significant effect on the gas boundary layer thickness (Aliseda et al. Reference Aliseda, Hopfinger, Lasheras, Kremer, Berchielli and Connolly2008). This equation is equivalent to that used in Varga et al. (Reference Varga, Lasheras and Hopfinger2003), where the parameter $b_g$ and the prefactor 2 are included in the constant. The relative velocity between the gas stream and the instability wave is defined as $u_r = u_g - u_l$. The constant $C_{\lambda }$ is determined by comparing (4.1) to measurements of the wavelength of the varicose instability, as detailed in Appendix A. For most practical cases, $u_l\ll u_g$, thus, $\lambda _{KH}$ and $u_r$ are relatively insensitive to the liquid flow rate; thus, this instability alone is insufficient to capture the effects of the changing liquid flow rate. Furthermore, while the dynamic of the liquid jet instability causes the liquid jet to be exposed directly to the air flow, it does not necessarily determine the diameter of the jet when it is atomized nor its relative velocity to the air stream, which are the governing parameters of the atomization.

To properly capture the effect of the liquid flow rate on the atomization, we describe a phenomenon that we call ‘liquid stream thinning’, where the low-speed liquid jet is accelerated to a higher speed by the air stream causing it to thin due to conservation of mass, as illustrated in figure 16. For a higher liquid flow rate, the jet will accelerate less, resulting in less thinning and a larger liquid jet diameter, ultimately leading to slightly larger breakup sizes. This phenomenon is visible in some of our images, although it is obscured by the instability and the breakup (for instance, figure 1); however, it is more clearly visible in the numerical results of Odier et al. (Reference Odier, Balarac and Corre2018) (see their figure 19), where the diameter of the jet decreases in the streamwise direction as its velocity increases. Although clearly present in their results, Odier et al. (Reference Odier, Balarac and Corre2018) did not comment on the effect as their work was concerned primarily on the nature of the flapping instability. As will be highlighted in future sections, the inclusion of such an effect offers a significant improvement in modelling the effect of the liquid flow rate on the breakup. Here, we consider a simple model for the liquid stream thinning, where it is assumed that the liquid jet accelerates from $u_l$ to the convective velocity, $u_c$, given by

(4.2)\begin{equation} u_c = \frac{\sqrt{\rho_l} u_l + \sqrt{\rho_g} u_g}{\sqrt{\rho_l} + \sqrt{\rho_g}}. \end{equation}

Figure 16. Illustration of the liquid stream thinning, where the liquid jet of initial diameter $d_l$ wicks across the face of the nozzle to $d_{g,i}$ and is accelerated by the surrounding gas flow, $u_g$, from $u_l$ to $u_c$, leading to the thinning of the liquid stream to $d_l'$.

The convective velocity is derived based on a balance of the dynamic pressures of the air and liquid flows on a wave moving in a Galilean frame of reference at $u_c$ (Dimotakis Reference Dimotakis1986). Therefore, by conservation of mass, the thinned liquid stream diameter, $d_l'$, is given by

(4.3)\begin{equation} d_l' = d_{g,i} \sqrt{\frac{u_l}{u_c}}. \end{equation}

The validity of this model component will be discussed further in § 4.5. We posit that the effect of the liquid flow rate on the atomization comes, primarily, from this liquid thinning phenomenon, where a higher $u_l$ will result in a slightly larger liquid jet diameter in the bulk sinuous waves than for a lower $u_l$. Another consequence of this description is that higher $u_g$ results in a higher degree of stream thinning, providing a slight amplification of the existing effects of increasing $u_g$.

Using both the physical and flow properties of the liquid and gas flows, as well as the description of the liquid jet's initial interaction with the gas jet (§ 4.1), the parameters required for defining the breakup of the jet based on the droplet breakup model can be determined; the liquid jet diameter, $d_l'$, Weber number, $We = \rho _g u_r d_l' / \sigma$, and time constant, $\tau = d_l' \sqrt {\rho _l / \rho _g}/ u_r$.

4.2. Breakup

Many of the existing phenomenological models for atomization of liquid jets use theories developed for the breakup of droplets to model the breakup of the waves in the liquid jet; most commonly, the Rayleigh–Taylor piercing model of Joseph et al. (Reference Joseph, Belanger and Beavers1999) (Varga et al. Reference Varga, Lasheras and Hopfinger2003; Aliseda et al. Reference Aliseda, Hopfinger, Lasheras, Kremer, Berchielli and Connolly2008). However, this model is limited in that it only predicts a single size that is usually compared with the $d_{32}$ of the spray and also does not describe any of the dynamics of the droplet that occur prior to breakup. Recent advances in aerodynamic droplet breakup modelling have yielded far more complete predictions of the dynamics of the droplet breakup. In the present work, the authors’ prior work on modelling aerodynamic droplet breakup (Jackiw & Ashgriz Reference Jackiw and Ashgriz2021, Reference Jackiw and Ashgriz2022) is applied to the breakup of the waves to provide a prediction of the resulting size distribution that results from the breakup of twin-fluid sprays. The advantages of this model are that it provides a more complete prediction of the droplet breakup, comparing well to both the droplet deformation and breakup, and that it provides a prediction of the child droplet size distribution based on the superposition of several breakup mechanisms. The model is also consistent with several droplet breakup morphologies, including sheet thinning; thus, its use in spray modelling extends this consistency to the morphologies of twin-fluid sprays, as discussed in the introduction (§ 1). In this section the aerodynamic droplet breakup model and its equations are summarized.

The aerodynamic droplet breakup model consists of modelling three phases of the droplet breakup, as follows:

1. (i) the droplet deformation, which defines the rim, bag and undeformed core (§ 4.2.1);

2. (ii) the breakup of the rim (§ 4.2.2) and the bag (§ 4.2.3) into 11 characteristic sizes;

3. (iii) the breakup of the undeformed core, which itself breaks as a droplet (§ 4.2.4).

After the breakup has been modelled, where 11 characteristic sizes are modelled for each breakup event (the first breakup as well as the breakup of the subsequent undeformed cores), the predicted characteristic sizes are interpreted as a size distribution or as statistics such as the $d_{32}$ (§ 4.3). For more details on (i), see Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021), with modifications in Jackiw & Ashgriz (Reference Jackiw and Ashgriz2022), and for more details on (ii), (iii) and the interpretation of the sizes as a distribution, see Jackiw & Ashgriz (Reference Jackiw and Ashgriz2022). Note that the parameters of the droplet breakup model introduced in this section are unchanged from the previous work on droplet breakup with the exception of the acceleration term for the undeformed core, which is discussed in § 4.2.4. Since this model contains many parameters, some of which are empirical based on underlying droplet breakup phenomena, it is important to note that they are not altered to fit to the present case of twin-fluid sprays and are based solely on droplet breakup.

4.2.1. Initial deformation and the formation of the rim and undeformed core

The first step in modelling the breakup of the effective droplet is to model its deformation to predict the dimensions of its deformed geometry. The geometries of interest are those of the deformed droplet at the initiation time (subscript $i$) and at the instant before breakup (subscript $f$). As the droplet deforms under the aerodynamic forces, liquid flows from the centre of the droplet face towards its equator, forming a disk at the droplet's windward face. Owing to a balance between the inertia of the radial flow and surface tension, a rim is formed along the periphery of the disk having a thickness at the initiation time of

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

where the cross-stream dimension of the rim ($R_i$) is given by the empirical relationship

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

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

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

Here $d_0$ denotes the diameter of the initial droplet, which can take the form of either $d_l'$ in the case of the initial breakup of the liquid jet or $d_c$ in the case of an undeformed core.

At low $We$, the entire droplet becomes the windward disk; however, at higher $We$, a portion of the droplet is left undeformed by the time the windward disk has formed. This remaining portion is called the ‘undeformed core’. The volumes of the windward disk and undeformed core are given by

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

and

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

respectively.

The undeformed core is the primary cause of the varying morphologies of breakup. Depending on the size and the relative velocity between the undeformed core and the air flow, the undeformed core may undergo further breakup. When the undeformed core is of negligible size, then the deformed droplet is comprised entirely of the disk structure ($V_c/V_0 \lessapprox 0.1$; note that 0.1 is used in the present work instead of 0 to make the code implementation practical) and the centre of the disk is blown into a single bag, constituting the bag breakup morphology. At moderate undeformed core sizes ($0.1 \lessapprox V_c/V_0 \lessapprox 0.5$), the conditions of the undeformed core are insufficient to cause its breakup; thus, the bag is constrained by the additional mass at the centre of the drop, which is drawn into a stamen, constituting the bag and stamen breakup morphology. At $V_c/V_0 \gtrapprox 0.5$ and $We<80$ (note that the other cases only occur at $We<80$), then the undeformed core will undergo its own breakup, constituting the multi-bag breakup morphology. Finally, at $We>80$, the surrounding airflow is strong enough to draw the rim downstream, forming a bag-like membrane between the undeformed core and the rim, constituting the sheet-thinning morphology.

After the initiation time, the region constrained between the rim and the undeformed core gets blown out into a bag. As the bag grows, it forces the rim to grow, thus, by conservation of mass, the rim is forced to thin throughout the bag's growth until it reaches its final dimension, given by

(4.9)\begin{equation} \frac{h_f}{d_0} = \frac{h_i}{d_0} \sqrt{\frac{2R_i}{d_0} \frac{d_0}{2R_f}}, \end{equation}

where $R_f$ is related to the morphology of the breakup and the size of the bag when it bursts, $\beta _f$. The relationship between $R_f$ and $\beta _f$ are determined for each morphology based on geometric considerations, and are summarized in table 2. Note that these geometric considerations include the bag-type and sheet-thinning morphologies, which were proposed in the introduction (§ 1) to be analogous to the membrane and fibre-type morphologies, respectively.

Table 2. Relationships between final rim diameter, $2R_f$, and bag size at burst, $\beta _f$, and their conditions for the four breakup morphologies. The acronyms B, BS, MB and ST denote bag breakup, bag and stamen breakup, multi-bag breakup and sheet-thinning morphologies, respectively.

The bag growth is found by a force balance at the tip of the bag, which results in

(4.10)$$\begin{gather} \frac{\beta\left(t^{*}\right)}{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\\ \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{gather}$$

where $t^{*} = t - t_i$ ($t_i\approx 0.9 \tau$), and the bag volume, $V_b/V_0$, is estimated as

(4.11)\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}

To find $\beta _f$, (4.10) is evaluated at the burst time, $t^{\ast }_b$, determined based on the instability of the accelerating bag tip as in Vledouts et al. (Reference Vledouts, Quinard, Vandenberghe and Villermaux2016) as

(4.12)\begin{equation} t^{{\ast}}_{b}=\dfrac{\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_b \dfrac{8 \tau}{\sqrt{3 We}} \sqrt{\dfrac{V_{b}}{V_{0}}} \dfrac{\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}

where $C_b=9.4$ is a constant relating to the breakup criterion of the bag. While this relationship has been shown by our previous work (Jackiw & Ashgriz Reference Jackiw and Ashgriz2021, Reference Jackiw and Ashgriz2022) to correlate well with the breakup time and size of the bags in aerodynamic droplet breakup, it is worth noting that other mechanisms may govern the rupture of the bag, such as impurities, Marangoni effects and film drainage, as discussed in Jackiw & Ashgriz (Reference Jackiw and Ashgriz2022) and Villermaux (Reference Villermaux2020).

4.2.2. Breakup of the rim

Since the rim of the deformed droplet is relatively long lived and has a changing geometry during its life, it undergoes multiple dynamics as it fragments. The breakup of the rim can be considered as having two distinct modes of breakup: the formation and breakup of large node droplets, and the breakup of the remaining rim between the nodes.

When the rim first forms, it is immediately subject to multiple possible instabilities; however, since the rim is constantly deforming through this period, the instabilities are nonlinear and lead to the formation of large nodes on the rim that are connected by the remaining rim segments. Once the distinct nodes are formed, with the remaining rim between them, they remain essentially constant in size through the breakup of the rim to form the node child drops. To account for the remaining rim, it is assumed that the nodes form from a (volume) fraction, $n_N$, of a wavelength segment of the initiated rim. From conservation of mass with such a segment and a spherical droplet, the size of the node drops are given by

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

where $\lambda$ is the wavelength of the instability that forms the node drops.

Both the Rayleigh–Plateau and Rayleigh–Taylor instability mechanisms give comparable results for the breakup of the rim. However, although the Rayleigh–Plateau mechanism seems phenomenologically more realistic and compares better with intermediate results for the number of nodes formed and the characteristic breakup time, the Rayleigh–Taylor instability gives a slightly better prediction of the node drop size, which is of primary importance in the present work. Zhao et al. (Reference Zhao, Liu, Li and Xu2010) derived the Rayleigh–Taylor wavelength for the rim of the deformed droplet as

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

where $C_D=1.2$ is the drag coefficient, and $2R_{max}$ is the maximum cross-sectional deformation at which the instability takes hold, which occurs slightly after the rim is initiated and is given by the following empirical correlation (Zhao et al. Reference Zhao, Liu, Li and Xu2010):

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

In (4.13), $n_N$ was determined empirically to have minimum, median and maximum values of $n_N=0.2, 0.4, 1$, respectively. For the purposes of the model, each value is taken such that the rim nodes provide three characteristic breakup sizes that capture the entire range of possible sizes. The relative-volume contribution of each of the node droplet characteristic sizes is determined as

(4.16)\begin{equation} w_N = \left( \frac{V_N}{V_0} \right) /3 \approx \left( 0.4 \frac{V_d}{V_0} \right) /3, \end{equation}

where the mean value of $V_N/V_d \approx 0.4$ is used for simplicity, and $V_d/V_0$ is given by (4.7). The factor $1/3$ accounts for the three different characteristic sizes that contribute to this mode of the breakup. Note that this factor was not included in the original analysis (Jackiw & Ashgriz Reference Jackiw and Ashgriz2022) due to a different treatment of the characteristic sizes when predicting the distribution, as will be described in § 4.3.

The remaining segments of the rim undergo breakup by the Rayleigh–Plateau instability; however, this instability manifests as two distinct mechanisms. In the conventional Rayleigh–Plateau instability, small perturbations on the rim, with a preferential wavelength set by the rim's minor diameter, grow until they fragment the rim, forming droplets of size (Ashgriz Reference Ashgriz2011)

(4.17)\begin{equation} \frac{d_{r}}{d_0} = 1.89 \frac{h_f}{d_0}. \end{equation}

An additional form of the Rayleigh–Plateau instability is manifested when the receding rim of the bag (discussed in § 4.2.3), which has a corrugation wavelength of $\lambda _{rr}$ (4.25, which will be discussed in § 4.2.3), collides with the remaining rim. The result is that the remaining rim receives a strong perturbation at the receding rim wavelength and will break at this wavelength instead of the most susceptible Rayleigh–Plateau wavelength when the collision is strong enough. Consequently, an additional characteristic droplet size is produced due to this collision mechanism that is given by

(4.18)\begin{equation} \frac{d_{r}}{d_0} = \left[ \frac{3}{2} \left( \frac{h_f}{d_0} \right)^2 \frac{\lambda_{rr}}{d_0} \right]^{1/3}. \end{equation}

Additionally, each of these mechanisms is accompanied by the formation of satellite droplets, whose sizes are given by Keshavarz et al. (Reference Keshavarz, Houze, Moore, Koerner and McKinley2020) as

(4.19)\begin{equation} d_s = d_r \frac{1}{\sqrt{2 + 3 Oh_{r}/\sqrt{2}}}, \end{equation}

where $Oh$ is the Ohnesorge number of the rim, given by $Oh_{r} = \mu _l / \sqrt {\rho _l h_f^3 \sigma }$. This relationship is based on an estimation of the thickness of the filament that is formed and stretched in the necking region between two wave crests in Rayleigh–Plateau breakup. Since $Oh_{r} \ll 1$ for the remaining rim, $d_s \approx d_r / \sqrt {2}$.

The relative-volume contribution of each of the remaining rim child droplet characteristic sizes is determined as

(4.20)\begin{equation} w_{r} = \left( \frac{V_r}{V_0} \right) /4, \end{equation}

where the rim volume, $V_r/V_0$, is given by

(4.21)\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}

As with the node drops, the factor $1/4$ in (4.20) accounts for the four characteristic sizes of the remaining rim's breakup: the Rayleigh–Plateau mechanism ($d_{r,RP}$, (4.17)), the collision mechanism ($d_{r,coll}$, (4.18)) and the accompanying satellite droplets for each ($d_{r, RP-s}$ and $d_{r, coll-s}$, (4.19)).

4.2.3. Breakup of the bag

The final and smallest characteristic sizes of the breakup are due to the breakup of the bag. When the bag perforates, it quickly recedes away from the perforation at the Taylor–Culick velocity (Culick Reference Culick1960),

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

where the minimum bag thickness is determined empirically as ${h}_{min} \approx 2.3\,\mathrm {\mu }$m. As it recedes, it follows the curvature of the bag, $\beta _f$, and experiences a centripetal acceleration,

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

and forms a rim governed by the universal rim thickness criterion of Wang et al. (Reference Wang, Dandekar, Bustos, Poulain and Bourouiba2018) ($Bo=\rho _l b_{rr}^2 a / \sigma = 1$). The receding rim thickness, $b_{rr}$, is found as

(4.24)\begin{equation} b_{rr} = \sqrt{\frac{\sigma}{\rho_l a_c}}. \end{equation}

The receding rim is then unstable to the Rayleigh–Plateau instability as

(4.25)\begin{equation} \lambda_{rr} = 4.5 b_{rr}, \end{equation}

leading to its fragmentation at a size of

(4.26)\begin{equation} \frac{d_{b,RP}}{d_0} = 1.89 \frac{b_{rr}}{d_0}, \end{equation}

with accompanying satellite droplets whose characteristic size are given by (4.19). Although this rim recession mechanism is clearly viewed in the breakup of the bag, the dynamics of the bag's breakup are the least well studied and understood of aerodynamic droplet breakup mechanisms. As such, two additional characteristic sizes are introduced to capture the range of sizes produced by the breakup of the bag. These additional characteristic sizes are the bag thickness, $h_{min}$, and the receding rim thickness, $b_{rr}$.

The relative-volume contribution of each of the remaining rim child droplet characteristic sizes is determined as

(4.27)\begin{equation} w_{b} = \left( \frac{V_b}{V_0} \right) /4 = \left( \frac{V_d}{V_0} - \frac{V_N}{V_0} - \frac{V_r}{V_0} \right) /4. \end{equation}

Once again, the factor $1/4$ in (4.27) accounts for the four characteristic sizes of the bag's breakup: the Rayleigh–Plateau mechanism on the receding rim ($d_{b,RP}$, (4.26)), its accompanying satellites ($d_{b,RP-s}$, (4.19)), the bag thickness (${h}_{min} \approx 2.3\,\mathrm {\mu }$m) and the receding rim thickness ($b_{rr}$, (4.24)).

4.2.4. Breakup of the undeformed core

After the first breakup of the effective droplet, there may remain an undeformed core, which itself may undergo further breakup (i.e. additional breakup events). Although the breakup model is still applicable, the input conditions for the model based on the undeformed core must be determined.

Within the subsequent breakup events, it is necessary to account for the volume weighting of the core in the weights of its characteristic sizes, i.e. the relative weights of each characteristic size of a core breakup event must be multiplied by the relative weight of the core breakup event itself, as the core breakup will constitute a smaller volume than the initial effective droplet. The relative-volume contribution of the undeformed core is determined as

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

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

(4.29)\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 (4.28). In the case of twin-fluid sprays, the gas speed is typically high enough that the change in speed of the undeformed core relative to the gas speed is negligible, thus, $u_r$ remains unchanged when modelling the breakup of the core. Using $u_r$ and the newly determined core droplet size, $d_c$, the core Weber number, $We_c$, and characteristic breakup time, $\tau _c$, can be determined, and the breakup model can be iterated until no further core breakup occurs ($We_c < 8.8$ or $V_c/V_0 < 0$).

Although the change in speed of the core droplets is not modelled, this does not mean that they break up immediately; in reality, they will essentially manifest as a breakup cascade as they break at different axial locations in the spray. It is postulated that this mechanism is responsible for the delay in transition of the bimodal distribution to a monomodal one in the axial direction described in § 3. Considering figure 7, the large mode would be attributed to yet unbroken undeformed cores, which, as they break, are transferred to the smaller sizes of the smaller mode. If the gas flow speed is sufficiently low, or the liquid flow speed sufficiently high, the core drops will be present further along the spray and may penetrate past the region where the air flow remains strong and will not break further, resulting in a bimodal distribution of droplet sizes. If the core drops were tracked in space by the model, or by a numerical scheme implementing the model, in addition to the falling speed of the expanding gas jet, the modality of the distribution could be predicted.

4.3. Prediction of the droplet size distribution

The models described in § 4.2 provide predictions of 11 unique characteristic sizes for each breakup event, which are summarized in table 3. Figure 17 shows how the set of characteristic sizes are constructed sequentially through the breakup of the undeformed cores by giving the size count of various breakup events where $j$ indexes the breakup event. Here $j=0$ is the first breakup event (panel a), $j=1$ is the breakup of the first undeformed core (panel b) up to the last breakup of an undeformed core at $j=18$ (panel c). The last breakup event leaves an undeformed core of size $d_c=8$ mm at $We_c = 6.5$, which is below the critical $We$ for breakup of 8.8 (Jackiw & Ashgriz Reference Jackiw and Ashgriz2021). Panel (d) gives the total size histogram of the entire breakup, consisting of a total of 209 droplets from 19 breakup events, while (e) gives the volume-weighted size histogram for the same. The input sizes and $We$ for each step are indicated on the figure.

Figure 17. Example construction of the predicted distribution for the 2850-70 nozzle at $G_l = 263$ and $G_g = 289$ kg m$^{-2}$ s$^{-1}$. (a–c) Size histograms of the 11 characteristic sizes for the breakup of the initial liquid jet ($j=0$) and the first ($j=1$) and last undeformed core ($j=18$), respectively. The conditions of each breakup event are indicated on the figure. (d) The total size and (e) volume-weighted histograms of the breakup. (f–h) The number and volume probability density distributions, and the volume-frequency histogram, respectively.

Table 3. List of characteristic breakup sizes.

The predicted characteristic sizes represent a simplification of the breakup dynamics to make the problem tractable, and as a result, their use alone will not provide a prediction of the distribution's shape. Therefore, the set of characteristic sizes requires interpretation through the figurative lens of a distribution density function to provide a prediction of the distribution, as shown by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2022) for aerodynamic droplet breakup.

While the authors’ prior model for droplet breakup (Jackiw & Ashgriz Reference Jackiw and Ashgriz2022) utilized a combination of Gamma distributions to predict a multimodal size distribution, in the present work, the sizes resulting from each of these modes is close enough that their contributions overlap to appear as a single distribution. Since the treatment of the undeformed core results in a dynamic somewhat similar to a cascade (since each core results from the prior breakup event), the use of a log-normal distribution is reasonably appropriate (Frisch & Sornette Reference Frisch and Sornette1997; Villermaux Reference Villermaux2007) (although, it is noted that there is continuing debate over what the ‘correct’ distribution function is in the representation of sprays). Furthermore, the data from the Malvern instrument is already at least somewhat inherently biased towards a log-normal distribution (or a combination of log-normal distributions), due to the way that the raw scattering data are corrected and interpreted (Canals et al. Reference Canals, Wagner, Browner and Hernandis1988; Wittner, Karbstein & Gaukel Reference Wittner, Karbstein and Gaukel2018; Sijs et al. Reference Sijs, Kooij, Holterman, Van De Zande and Bonn2021) (note that this seems to be the case even in the software's ‘model independent’ mode). It is therefore appropriate to model the distribution as log normal for the purposes of comparing to such data. As such, it can be assumed that the number-density distribution, $p_n$, of the spray follows a log-normal probability density distribution, as

(4.30)\begin{equation} p_n(x=d) = \frac{1}{x s \sqrt{2 {\rm \pi}}} \exp \left( - \frac{\left(\ln(x)- \bar{x} \right)^2}{2 s^2} \right), \end{equation}

where $\bar {x}$ and $s$ are the distribution's mean and standard deviation, classically given by

(4.31a,b)\begin{equation} \bar{x} = \ln \left( \frac{\bar{x}_X^2}{\sqrt{\bar{x}_X^2 + s_X^2}} \right), \quad s^2 = \ln \left( 1 + \frac{s_X^2}{\bar{x}_X^2} \right), \end{equation}

where $\bar {x}_X$ and $s_X$ are the weighted mean and standard deviation of the set of predicted breakup sizes weighted by their respective modes, given by

(4.32a,b)\begin{equation} \bar{x}_X=\dfrac{\sum_{i=1}^{\check{n}} w_{i} x_{i}}{\sum_{i=1}^{\check{n}} w_{i}}, \quad s_X = \sqrt{\frac{\sum_{i=1}^{\check{n}} w_{i}\left(x_{i}-\bar{x}_{X}\right)^{2}}{\dfrac{(\check{n}-1)}{\check{n}} \sum_{i=1}^{\check{n}} w_{i}}}, \end{equation}

$x_i$ are the characteristic sizes of each mode ($x=d$) and $\check {n}$ is the total number of characteristic sizes. An example of the resulting $p_n$ is shown in figure 17(f).

The $p_n$ of a droplet size distribution can often be misleading, as small droplets will dominate the $p_n$ due to their large numbers even if they contain a negligible fraction of the spray's volume. It is therefore prudent to compare the volume-weighted probability density of the spray, $p_v$, rather than the $p_n$, as this allows for a more realistic view of how the volume of the spray is distributed across the spray sizes, which is more useful in most applications. From the $p_n$, the $p_v$ is obtained by weighting the function by $x^3$ and renormalizing the distribution, as

(4.33)\begin{equation} p_v(x) = \frac{x^3 p_n(x)}{\int_{0}^{\infty} x^3 p_n(x)\,{{\rm d}\kern0.7pt x} } = p_n \left(\frac{x^3}{\exp \left( 3 \bar{x} + 4.5 s^2 \right)} \right). \end{equation}

An example of the resulting $p_v$ is shown in figure 17(g).

Finally, to compare directly to the results of the Malvern Spraytec, the $p_v$ must be converted to the volume frequency, $f_v$, which requires the distribution to be integrated over discrete bins to obtain the frequency distribution from the probability density distribution, as

(4.34)\begin{equation} f_v \left(x_{i,1} \rightarrow x_{i,2} \right) = \int_{x_{i,1}}^{x_{i,2}} p_v(x) \, {{\rm d}\kern0.7pt x} \; 100\,\%, \end{equation}

where $x_{i,1}$ and $x_{i,2}$ are the left and right edges of the $i$th bin, respectively. For direct comparison of the Malvern Spraytec measurements, the same bins as the Malvern results are used in the conversion. An example of the resulting $f_v$ is shown in figure 17(h). Strictly speaking, (4.34) gives the total probability of each bin, not the frequency. The frequency represents an actual result (e.g. the measurement of the instrument), whereas the probability is the theoretical analogue. For simplicity and consistency in the plotting, they are both labelled as $f_v$ in this text.

It should be pointed out that even though the predicted characteristic sizes do not exceed 25 $\mathrm {\mu }$m in the case presented in figure 17 (see panels d,e), which is also reflected in the $p_n$ (panel f), while the predicted $p_v$, and therefore $f_v$, show sizes up to $100\,\mathrm {\mu }$m (see panels g,h). The continuous number distribution used, whose parameters were determined from the predicted characteristic sizes, theoretically extends to infinite $d$ and has small but non-zero values above $25\,\mathrm {\mu }$m that are amplified in the conversion to the volume-weighted distribution. It is therefore important to emphasize that the predicted sizes are simply characteristic of the breakup, which is why they are sufficient for estimating the distribution parameters, but not necessarily for directly representing the distribution.

The model predictions of the droplet size distribution are compared with the measured distributions of the spray from the 2850-70 nozzle at varying $G_g$ and at low and high $G_l$ in figure 18. To clarify comparison between the predicted and measured $f_v$, the predicted $f_v$ is presented as a solid line where dot markers are placed at the bin centres to highlight the bin heights. At low $G_l$, the prediction matches the measured distribution well for all cases with the exception of at low $G_g$ (see figure 18a) where the atomization is poor. Since the experimental distribution is bimodal at this condition, it has an excess of sizes ${>}100\, {\rm \mu}$m that is not captured by the present model; however, the main peak, which results from the aerodynamic breakup, is predicted well. Note that since the predicted distribution is slightly narrower than the experimental distributions, the predicted frequency near the main peak is overall higher. At high $G_l$, the model underpredicts the sizes of the distribution, which indicates that the effect of the liquid flow rate is not yet captured correctly. This was alluded to in § 4.1 when the liquid stream thinning effect was introduced, and will be discussed further in § 4.5. However, it is worth noting that the presented high $G_l$ case is substantially higher than the low $G_l$ case compared with the range of $G_g$, highlighting that the effect of the liquid flow rate on the atomization is relatively small. Similar to the low $G_l$ cases, when the size distribution is bimodal, the model is closer to the sizes that result from aerodynamic breakup and the large mode that results from unbroken portions of the liquid jet is not predicted. Similar results are found for the larger 60100-120 nozzle for varying $G_g$ and for low and high gas $G_l$, as shown in figure 19, demonstrating that the model performs consistently across both nozzle scales.

Figure 18. Comparison of predicted and measured droplet size distribution for the 2850-70 nozzle at low (a–d) and high (e–h) $G_l$ and increasing $G_g$ at $x=80$ mm. Flow conditions are indicated on the figure. Results are shown for (a) $G_g = 143$, $G_l = 263\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (b) $G_g = 210$, $G_l = 263\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (c) $G_g = 289$, $G_l = 263\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (d) $G_g = 535$, $G_l = 263\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (e) $G_g = 143$, $G_l = 1315\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (f) $G_g = 210$, $G_l = 1315\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (g) $G_g = 289$, $G_l = 1315\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; and (h) $G_g = 535$, $G_l = 1315\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$.

Figure 19. Comparison of predicted and measured droplet size distribution for the 60100-120 nozzle at $G_l=164$ kg m$^{-2}$ s$^{-1}$ and increasing $G_g$ from (a–d). Flow conditions are indicated on the figure. Results are shown for (a) $G_g = 140$, $G_l = 164\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (b) $G_g = 200$, $G_l = 164\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (c) $G_g = 328$, $G_l = 164\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (d) $G_g = 527$, $G_l = 164\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (e) $G_g = 140$, $G_l = 493\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (f) $G_g = 200$, $G_l = 493\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (g) $G_g = 328$, $G_l = 493\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (h) $G_g = 527$, $G_l = 493\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$.

4.4. Prediction of the mean diameter and comparison to other models

Alternatively, for a simpler representation of the model, the $d_{32}$ can be calculated from the estimated number probability density as

(4.35)\begin{equation} d_{32} = \frac{\int_0^{\infty} x^3 p_n(x)\, {{\rm d}\kern0.7pt x}}{\int_0^{\infty} x^2 p_n(x) \,{{\rm d}\kern0.7pt x}}. \end{equation}

Quantifying the model performance based on $d_{32}$ has two main advantages. Firstly, it provides a straightforward method to test the model's performance over a wide range of parameters. Secondly, since most prior models (that do not require the assumption of an unknown level of variability) offer a prediction only for $d_{32}$, this is the only way to effectively compare the present model to existing ones. In this section the present model, using (4.35) to predict $d_{32}$, will be compared with the models of Varga et al. (Reference Varga, Lasheras and Hopfinger2003) (A1) and Aliseda et al. (Reference Aliseda, Hopfinger, Lasheras, Kremer, Berchielli and Connolly2008) (A2). The determination of the empirical coefficients for these models is described in Appendix A.

The present model along with those of Varga et al. (Reference Varga, Lasheras and Hopfinger2003) and Aliseda et al. (Reference Aliseda, Hopfinger, Lasheras, Kremer, Berchielli and Connolly2008) are compared with the experimentally measured $d_{32}$ of the 2850-70 nozzle sprays for varying $G_g$ and for low $G_l$ in figure 20(a). At these flow conditions, the present model and that of Varga et al. (Reference Varga, Lasheras and Hopfinger2003) provide nearly identical results, while the model of Aliseda et al. (Reference Aliseda, Hopfinger, Lasheras, Kremer, Berchielli and Connolly2008), although close, does not match the trend at high $G_g$. The prediction of the present model alone is compared with the experimentally measured $d_{32}$ in figure 20(b), showing excellent agreement across the range. Figure 20(b) also highlights the cases that exhibit a bimodal distribution (open markers).

Figure 20. (a) Comparison of the present model along with those of Varga et al. (Reference Varga, Lasheras and Hopfinger2003) and Aliseda et al. (Reference Aliseda, Hopfinger, Lasheras, Kremer, Berchielli and Connolly2008), and (b) comparison of the present model alone to the measured $d_{32}$ of the 2850-70 nozzle spray for varying $G_g$ at low $G_l$ and $x=80$ mm. Same for (c,d) at high $G_l$. Flow conditions are indicated on the figure. Results are shown for (a,b) $G_l =263\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$, (c,d) $G_l =1315\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$.

The same is shown for high $G_l$ in figure 20(c,d). Note that at very high $G_l$ and low $G_g$, the strong large mode of the bimodal distribution causes a significant increase in the experimental $d_{32}$, which is evidenced by the rightmost datum in figure 20(d) that deviates from the overall trend of the data. At these conditions, both the present model and that of Varga et al. (Reference Varga, Lasheras and Hopfinger2003) underpredict the measurement, indicating that the effect of the liquid flow rate is not properly captured in either. Meanwhile, the model of Aliseda et al. (Reference Aliseda, Hopfinger, Lasheras, Kremer, Berchielli and Connolly2008), although still not capturing the trend with $G_g$ properly, does exhibit a more accurate sensitivity to $G_l$. These variations highlight the differences in how each model incorporates the effect of the liquid flow rate. All of the models include the effect of $u_l$ on the atomization in the relative speed between the gas and liquid flows, $u_r = u_g - u_l$; however, as mentioned in § 4.1, since $u_l\ll u_g$ for most practical cases, including the effect of $u_l$ on this term alone will not be sufficient to capture the effects of the liquid flow rate on the atomization.

In an attempt to capture the liquid flow effect, Aliseda et al. (Reference Aliseda, Hopfinger, Lasheras, Kremer, Berchielli and Connolly2008) add the familiar empirical term, $(1+m_r)$, which arises from energy considerations, to their model. While this term has been shown to be relevant in many twin-fluid atomization cases, its addition to the model in this way is naive of the actual phenomenology that gives rise to it. To provide a phenomenological explanation of the effect of the liquid flow rate on the atomization, the present model proposes the liquid stream thinning effect, where the advection of the liquid stream affects its diameter and in turn the characteristic size that is exposed to the gas flow to undergo breakup (see § 4.1). Since the present model still underpredicts the sensitivity of the atomization to $G_l$, the proposition of the liquid stream thinning effect appears questionable; however, it may merely not yet be properly modelled.

Similar results are found for the larger 60100-120 nozzle, as shown in figure 21. At the larger nozzle scale, the present model captures the effect of $G_l$ on the atomization slightly better than that of Varga et al. (Reference Varga, Lasheras and Hopfinger2003). Note that at the high $G_l$ case (figure 21cd), all but the highest two $G_g$ cases exhibit bimodal size distributions.

Figure 21. (a) Comparison of the present model along with those of Varga et al. (Reference Varga, Lasheras and Hopfinger2003) and Aliseda et al. (Reference Aliseda, Hopfinger, Lasheras, Kremer, Berchielli and Connolly2008), and (b) comparison of the present model alone to the measured $d_{32}$ of the 60100-120 nozzle spray for varying $G_g$ at low $G_l$ and $x=80$ mm. Same for (c) and (d) at high $G_l$. Flow conditions are indicated on the figure. Results are shown for (a,b) $G_l =164\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$, and (c,d) $G_l =493\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$.

4.5. The effect of liquid flow rate

As demonstrated in the previous section, the effect of increasing the liquid flow rate has not been completely captured by the present model. The liquid stream thinning model (4.3) assumes that the entirety of the jet is accelerated to the convective velocity, $u_c$ (4.2), by the time it undergoes breakup; however, there are two main issues with this assumption. Firstly, as described in the introduction of $u_c$ in § 4.1, the convective velocity is based on the balance of the stagnation pressures of the air and liquid flows of a surface wave in a Galilean frame of reference (Dimotakis Reference Dimotakis1986). This description essentially assumes a complete transfer of momentum to the liquid wave from the air stream. However, the nature of the described liquid stream thinning effect is that the liquid jet is dragged tangentially by the air stream to a higher speed, in which case the momentum transfer is not complete. Secondly, the described liquid stream thinning only occurs when the liquid jet is advecting differentially along its axis. Once a portion of the liquid jet is advecting as one, no liquid stream thinning will occur. Such a case occurs when the jet's instability transitions into a flapping or helical mode where the jet is accelerated nearly perpendicular to its axis; thus, the entire wave is advected as a whole and no liquid stream thinning occurs. Since the liquid stream thinning only occurs up to this point, it will not occur all the way to the convective velocity, $u_c$, as modelled.

A more realistic description of the liquid stream thinning, therefore, is the thinning that occurs due to the advection of the liquid jet by the air stream until it undergoes the flapping instability. While only subtly different from the original description, this modified description crucially includes the provision that the liquid jet will not have necessarily fully accelerated to some terminal velocity by the time the thinning ceases to occur. To provide a model for this modified depiction of the liquid stream thinning, we assume a different scenario for the determination of the liquid jet speed at the time of interest, wherein momentum is transferred incompletely to the liquid jet from the airflow. Considering the case where only some fraction, $C_s$, of the air jet's momentum is transferred to the liquid jet, we obtain the accelerated liquid velocity as

(4.36)\begin{equation} u_{l}' = \sqrt{u_l^2 + C_s \left( \frac{\rho_g}{\rho_l} \right) u_g^2}, \end{equation}

where $C_s\approx 0.1$ gives good agreement for the overall model and is physically reasonable (i.e. approximately 10 % of the gas flow momentum is transferred to the liquid jet by the time the liquid stream develops bulk waves that advect as a whole).

The results of the modified model, now using (4.36) instead of (4.2) in (4.3), are compared with the present experiments for the size distribution in figure 22 and the $d_{32}$ in figure 23 for the 2850-70 nozzle. In both figures it is evident that the modified model provides an improved prediction for the higher $G_l$ cases (compare to figures 18 and 20 for the size distribution and $d_{32}$, respectively). However, this model is mainly an empiricism to account for the dynamic of the jet prior to its sinuous instability, which has been neither accounted for nor studied previously. Further study of the early interaction between the liquid and air flows will be required to properly understand and model the effects of the liquid flow rate. Another factor that may further complicate the thinning of the liquid jet is the stretching of the jet as the amplitude of the flapping instability increases; however, this is omitted from the present analysis.

Figure 22. Comparison of predicted and measured droplet size distribution for the 2850-70 nozzle at low (a–d) and high (e–h) $G_l$ and increasing $G_g$ at $x=80$ mm. Flow conditions are indicated on the figure. The prediction uses $u'_l$ (4.36) in the prediction of the liquid stream thinning. Results are shown for (a) $G_g = 143$, $G_l = 263\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (b) $G_g = 210$, $G_l = 263\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (c) $G_g = 289$, $G_l = 263\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (d) $G_g = 535$, $G_l = 263\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (e) $G_g = 143$, $G_l = 1315\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (f) $G_g = 210$, $G_l = 1315\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (g) $G_g = 289$, $G_l = 1315\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$; (h) $G_g = 535$, $G_l = 1315\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$.

Figure 23. (a) Comparison of the present model along with those of Varga et al. (Reference Varga, Lasheras and Hopfinger2003) and Aliseda et al. (Reference Aliseda, Hopfinger, Lasheras, Kremer, Berchielli and Connolly2008), and (b) comparison of the present model alone to the measured $d_{32}$ of the 2850-70 nozzle spray for varying $G_g$ at low $G_l$ and $x=80$ mm. Same for (c,d) at high $G_l$. Flow conditions are indicated on the figure. Results are shown for (a,b) $G_l =263\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$, (c,d) $G_l =1315\,{\rm kg}\,{\rm m}^{-2}\,{\rm s}^{-1}$.

4.6. Effects of phase properties on the distribution

So far, the model has been validated over a practical range of gas and liquid flow rates and nozzle sizes for air-assisted water sprays only. While the effects of the gas and liquid properties on the atomization is of interest in many applications, comprehensive reporting of the size distribution of sprays for other fluids is limited, impeding the validation of the present model over such ranges. To assess the theoretical effects and sensitivity of varying the phase properties on the atomization, the present model is computed over a wide and practical range of fluid properties, as presented in figure 24.

Figure 24. Model results over several orders of magnitude of (a) $\rho _g$, (b) $\rho _l$ and (c) $\sigma$. In (a), $u_g = 100$ m s$^{-1}$ and $\rho _g=1.24$ kg m$^{-3}$, while in (b) and (c), $u_g = 200$ m s$^{-1}$ and $\rho _g=1.36$ kg m$^{-3}$. In all cases, $u_l = 0.3$ m s$^{-1}$. Unless otherwise specified, the liquid properties used are $\rho _l = 1000$ kg m$^{-3}$ and $\sigma =0.0729$ N m$^{-1}$, which are the nominal properties for the water base case at similar conditions.

In figure 24(a), $\rho _g$ is varied over a range of $\rho _g=1\unicode{x2013}300$ kg m$^{-3}$. Gas densities $\rho _g$ of the order O(1–10) is typical for many industrial manufacturing applications, such as spray drying, while O(100) is typical for combustion applications, such as reciprocal and turbine engines. The size distribution is found to decrease with increasing $\rho _g$, which appears throughout many stages of the model and in all cases causes the reduction in sizes. The most dominant effect of $\rho _g$ is to increase the radial deformation rate (4.6), leading to the decrease of the initial rim thickness (4.4), which in turn results in lower sizes for the main breakup mechanisms (i.e. the nodes and the remaining rim).

In figure 24(b), $\rho _l$ is varied over a range of $\rho _l = 500\unicode{x2013}20\,000$ kg m$^{-3}$, where the size distribution is shown to increase with increasing $\rho _l$; however, the prediction is not nearly as sensitive to $\rho _l$ as it was shown to be to $\rho _g$. Liquid densities of the order O(500–5000) are typical for most applications, while high liquid densities of the order O(10 000) are applicable for the atomization of liquid metals (note, however, that the secondary effects that would typically accompany such fluids, such as high liquid viscosity and chemical reaction rate, are not modelled). While increasing $\rho _l$ slows the deformation rate (4.6), the effect on the prediction of the rim thickness cancels due to the rim's inertia (4.4). While $\rho _l$ appears in the derivation of terminal bag dimensions (4.10) and serves to decrease the thickness of the rim through extending the lifetime of the bag, the effect is relatively small. Additionally, $\rho _l$ cancels in the prediction of $b_{rr}$ (4.24) and so does not significantly affect the breakup of the rim due to the collision of the receding rim. Furthermore, $\rho _l$ appears in the satellite droplet size prediction of Keshavarz et al. (Reference Keshavarz, Houze, Moore, Koerner and McKinley2020) (4.19) and results in bigger satellites; however, for low viscosity cases, the $Oh$ is low and the effect is small. Indeed, the main effect of $\rho _l$ on the atomization is in the liquid stream thinning, where higher liquid density results in less thinning due to a lower $u_c$, thus, $d_l'$ is slightly bigger, resulting in a slightly larger distribution.

Following the discussion of the different effects of varying $\rho _g$ and $\rho _l$, it is worth noting that most works that consider varying $\rho _l$ and/or $\rho _g$, in particular those on droplet breakup (see, for example, Jain et al. Reference Jain, Tyagi, Prakash, Ravikrishna and Tomar2019; Marcotte & Zaleski Reference Marcotte and Zaleski2019; Joshi & Anand Reference Joshi and Anand2022), frequently only report their trends in terms of the density ratio, $\rho _l/\rho _g$. In most of these works, only $rho_l$ is varied, which may lead to the potentially erroneous conclusion that decreasing $\rho _g$ should have the same effect as increasing $\rho _l$, so long as the ratio is the same. However, the present results suggest that the atomization has a different sensitivity to each. As a result, reporting such trends in terms of the density ratio only may lead to conflicting conclusions, as studies that only vary $\rho _l$ will report low sensitivity to the density ratio, while studies that vary $\rho _g$ will report high sensitivity. While the density ratio may still play a role in the atomization (for instance, in the flapping instability, as reported in Matas, Delon & Cartellier Reference Matas, Delon and Cartellier2018), it does not uniquely determine the effects of the phase densities, as their absolute values are also important.

Figure 24(c) shows the effect of $\sigma$ on the distribution for a range of $\sigma = 0.01\unicode{x2013}1$ N m$^{-1}$. Here $\sigma$ on the order of O(0.01) is applicable for solvents (e.g. for acetone, $\sigma = 0.0252$ N m$^{-1}$), while the order of O(0.1–1) is applicable for liquid metals. Surface tension $\sigma$ appears throughout the model and typically serves to counteract the deformation forces. As a result, increasing $\sigma$ results in an increase in the predicted size distribution, with a sensitivity to $\rho _g$.