1. Introduction
Drop fragmentation, also known as secondary fragmentation, is the process through which a drop breaks up under the action of external aerodynamic forces induced by the ambient flow. These forces originate due to a velocity deficit between the drop and the ambient medium. There are two fundamental ways a drop might experience a velocity deficit: a uniform ambient flow impacts a stationary drop in a gravity-free environment, called ‘impulsive acceleration’ (Han & Tryggvason Reference Han and Tryggvason2001); or an initially stationary drop accelerates under the action of a constant body force, while experiencing aerodynamic forces, called ‘free fall’ (Jalaal & Mehravaran Reference Jalaal and Mehravaran2012). For both cases, a liquid drop experiences aerodynamic forces that cause the drop to deform, which may lead to its fragmentation at a Weber number
$\mathit{We}_0$
higher than the critical value
$\mathit{We}_{cr}$
(Hinze Reference Hinze1949, Reference Hinze1955). During free fall, the drop starts with zero aerodynamic forces that gradually increase to a maximum, at either its terminal or its breakup velocity. On the other hand, an impulsively accelerated drop starts its deformation process with the largest velocity deficit, and the corresponding large aerodynamic stresses acting on it. These stresses gradually reduce as the drop decelerates with respect to the ambient flow. It should be noted that as the drop decelerates, it also simultaneously deforms causing an increase in its frontal area, which can in turn increase surface shear stresses, given the velocity deficit is still substantial.
Most applications such as the internal combustion engine, spray painting, etc. involve a secondary fragmentation due to impulsive acceleration. Among impulsive acceleration cases, there can be different experimental systems such as a drop introduced into a uniform cross-flow or a drop exposed to a shockwave in a wind tunnel. The time scales of an impulsive drop breakup process are relatively small, resulting in both the aforementioned experimental set-up predicting a similar critical Weber number for secondary fragmentation (Hsiang & Faeth Reference Hsiang and Faeth1992, Reference Hsiang and Faeth1995). Most of the experimental studies conducted on secondary fragmentation (Pruppacher & Beard Reference Pruppacher and Beard1970; Krzeczkowski Reference Krzeczkowski1980; Wierzba Reference Wierzba1990; Hsiang & Faeth Reference Hsiang and Faeth1992; Gelfand Reference Gelfand1996; Theofanous, Li & Dinh Reference Theofanous, Li and Dinh2004; Szakáll et al. Reference Szakáll, Diehl and Mitra2009; Kulkarni & Sojka Reference Kulkarni and Sojka2014; Jain et al. Reference Jain, Prakash, Tomar and Ravikrishna2015) have focused on impulsive acceleration, and those results are summarised in figure 1.
(a--e) Types of drop breakup morphologies (Guildenbecher, López-Rivera & Sojka Reference Guildenbecher, López-Rivera and Sojka2009; Theofanous Reference Theofanous2011) observed in experiments in order of increasing threshold Weber numbers. Panel (f) plots all the experimental data on threshold Weber numbers required to produce different fragmentation morphologies, based on a similar plot in Hsiang & Faeth (Reference Hsiang and Faeth1995) and data from Krzeczkowski (Reference Krzeczkowski1980), Pilch & Erdman (Reference Pilch and Erdman1987), Wierzba (Reference Wierzba1990), Dai & Faeth (Reference Dai and Faeth2001), Han & Tryggvason (Reference Han and Tryggvason2001), Kulkarni & Sojka (Reference Kulkarni and Sojka2014), Jain et al. (Reference Jain, Tyagi, Prakash, Ravikrishna and Tomar2019), Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021). Backward bag breakup, as shown is panel (b) in the red box and the red curve in panel (f), has been the predominantly observed critical non-vibrational fragmentation morphology.
After starting from an initial spherical shape, drops start the deformation process with flattening of the downstream face under the influence of pressure forces (Villermaux & Bossa Reference Villermaux and Bossa2009; Jain et al. Reference Jain, Tyagi, Prakash, Ravikrishna and Tomar2019; Jackiw & Ashgriz Reference Jackiw and Ashgriz2021). This is followed by the formation of a pancake of one of the following two types: (i) a flat disk-like structure with both upstream and downstream faces showing an increase in radius of curvatures (henceforth called ‘flat pancake’); or (ii) a pancake with a concave-shaped downstream surface, corresponding with minimal change in curvature of the upstream surface (henceforth called ‘forward pancake’) (Han & Tryggvason Reference Han and Tryggvason2001). These differences in pancake shapes have been observed for differences in physical and flow parameters such as density ratio
$\rho$
and initial Reynolds number
$\mathit{Re}_0$
(or outside Ohnesorge number
$\mathit{Oh}_o$
). However, the exact physical mechanism that causes this difference in pancake morphology has not yet been explored in the literature. Beyond the formation of a pancake, the pancake deforms further and starts to form a toroidal periphery (rim), which then leads to further deformation and possibly even breakup through different morphologies (discussed in the next paragraph). This stage, which marks the completion of pancake formation and the start of a visible peripheral rim, can be temporally indicated by a non-dimensional time
$t^* = t/\tau \approx 1$
. Here,
$\tau = D\sqrt {\rho }/V_0$
is the drop deformation time scale (Rimbert et al. Reference Rimbert, Castrillon Escobar, Meignen, Hadj-Achour and Gradeck2020), where
$V_0$
is the uniform initial velocity of the ambient medium relative to the drop,
$D$
is the drop’s volume-averaged diameter and
$t$
represents the elapsed dimensional time during the deformation process. This time scale is the same as the dimensionless time for Rayleigh–Taylor or Kelvin–Helmholtz instabilities specified by Pilch & Erdman (Reference Pilch and Erdman1987). Here
$\tau$
includes the effect of
$\rho$
, thus making
$t^*$
a useful temporal scale when comparing cases with different density ratios.
Following the formation of a pancake after
$t^*\gt 1$
, the drop may further deform and ultimately breakup through one of the following morphologies (as illustrated in figure 1
a–e): (i) the vibrational mode where the drop oscillates about a maximum deformation state, and does not show consistent breakup (Hsiang & Faeth Reference Hsiang and Faeth1992; Rimbert et al. Reference Rimbert, Castrillon Escobar, Meignen, Hadj-Achour and Gradeck2020); (ii) the simple bag breakup that involves the formation of a toroidal rim and the inflation of a thin film (bag) in between, which ultimately ruptures due to Rayleigh–Plateau instabilities (Kulkarni & Sojka Reference Kulkarni and Sojka2014; Jackiw & Ashgriz Reference Jackiw and Ashgriz2021); (iii) a bag breakup with morphological features in addition to a bag, such as a stamen/plume (Hsiang & Faeth Reference Hsiang and Faeth1995; Jain et al. Reference Jain, Prakash, Tomar and Ravikrishna2015) or multiple bags (Cao et al. Reference Cao, Sun, Li, Liu and Yu2007; Jackiw & Ashgriz Reference Jackiw and Ashgriz2021); (iv) sheet-thinning breakup where thin sheets and ligaments are removed from the periphery of the pancake and are blown downstream relative to the drop core due to their low local inertia, ultimately breaking up due to instabilities (Khosla & Smith Reference Khosla and Smith2006; Guildenbecher et al. Reference Guildenbecher, López-Rivera and Sojka2009); and (v) catastrophic breakup where unstably growing surface waves pierce through the entire pancake and cause it to catastrophically disintegrate (Theofanous Reference Theofanous2011).
For liquid drops in air under standard atmospheric conditions,
$\rho \gt 500$
and
$0.0005\lt \mathit{Oh}_o\lt 0.005$
. Consequently, in the context of density ratio, most existing experimental works are concentrated within
$\rho \gt 500$
. For all such experimental works, the critical drop breakup morphology for
$\mathit{Oh}_d\lt 0.1$
has been observed to always be a simple bag breakup, with minimal dependence of
$\mathit{We}_{cr}$
on density ratio. Through the advent of petascale/exascale computing in the last two decades, both axisymmetric and three-dimensional (3-D) direct numerical simulations of large
$\rho$
cases has become possible, thus allowing computational exploration of the entire density-ratio space for drops, i.e.
$\rho \gt 1$
. As a result, substantial research has been conducted in the high
$\rho$
space (Theofanous Reference Theofanous2011; Tavangar, Hashemabadi & Saberimoghadam Reference Tavangar, Hashemabadi and Saberimoghadam2015; Strotos et al. Reference Strotos, Malgarinos, Nikolopoulos and Gavaises2016; Yang et al. Reference Yang, Jia, Che, Sun and Wang2017; Guan et al. Reference Guan, Liu, Wen and Shen2018; Dorschner et al. Reference Dorschner, Biasiori-Poulanges, Schmidmayer, El-Rabii and Colonius2020; Ling & Mahmood Reference Ling and Mahmood2023; Tang, Adcock & Mostert Reference Tang, Adcock and Mostert2023), where threshold breakup morphologies similar to the experiments have been observed. On the other hand, computational works exploring low density-ratio cases (
$\rho \lt 100$
) tell a different story. Han & Tryggvason (Reference Han and Tryggvason2001) was one of the earliest computational works to explore the role of density ratio in the impulsive acceleration of a drop for various ambient and drop viscosities. The work focused on
$\rho =1.15$
and
$10$
, and found
$\rho$
to significantly influence the threshold Weber numbers and the corresponding fragmentation morphology (e.g. forward pancake and bag formation). Jain et al. (Reference Jain, Tyagi, Prakash, Ravikrishna and Tomar2019) explored the effect of
$\rho$
on drop deformation for a specific
$\mathit{Re}_0$
and viscosity ratio for a range of
$\mathit{We}_0$
from
$20$
to
$100$
, and observed the immense impact
$\rho$
has on the orientation of bags and pancakes, the drop velocities and the total observed deformations. Similar conclusions were reached by Marcotte & Zaleski (Reference Marcotte and Zaleski2019), where a higher threshold Weber number for both fragmentation and the transition from bursting to stripping was observed for the low density-ratio cases.
By the 1990s, experimental and theoretical studies (Karam & Bellinger Reference Karam and Bellinger1968; Krzeczkowski Reference Krzeczkowski1980; Pilch & Erdman Reference Pilch and Erdman1987) had established the important role of the drop Ohnesorge number
$\mathit{Oh}_d$
on
$\mathit{We}_{cr}$
. Hsiang & Faeth’s review paper in 1995 (Hsiang & Faeth Reference Hsiang and Faeth1995) significantly advanced this understanding. They aggregated all available experimental data from existing works, including their own experiments, into
$\mathit{We}_{cr}$
vs
$\mathit{Oh}_d$
plots. Their findings showed that the threshold
$\mathit{We}_0$
for the onset of all types of breakup morphologies (i.e. simple backward bag and other higher
$\mathit{We}_0$
breakup morphologies) follows a similar trend with respect to
$\mathit{Oh}_d$
(see figure 1 of Hsiang & Faeth (Reference Hsiang and Faeth1995) or figure 1
f), with the threshold
$\mathit{We}_0$
almost independent with respect to
$\mathit{Oh}_d$
for
$\mathit{Oh}_d\lt 0.1$
, and then increasing rapidly for
$\mathit{Oh}_d\gt 0.1$
. Furthermore, the critical breakup morphology (for the onset of breakup) for all the explored works were found to be simple bag breakups. It is crucial to emphasise that all the findings presented in Hsiang & Faeth (Reference Hsiang and Faeth1995) were derived from experiments in the high density-ratio regime, i.e.
$\rho \gt 500$
. This focus on high density ratios is also reflected in the majority of computational studies investigating the effect of
$\mathit{Oh}_d$
on drop breakup. For instance, Yang et al. (Reference Yang, Jia, Che, Sun and Wang2017) observed an increase in the transition Weber number from squeezing to bag breakup for increasing
$\mathit{Oh}_d$
for drops with
$\rho =800$
. Similarly, Jain et al. (Reference Jain, Tyagi, Prakash, Ravikrishna and Tomar2019) explored the effect of viscosity ratio (and hence
$\mathit{Oh}_d$
) on breakup morphologies through simulations of drops of
$\rho =1000$
and two different viscosity ratios. They noted a decrease in
$\mathit{We}_0$
with decreasing
$\mathit{Oh}_d$
and the emergence of a plume at the upstream pole for lower viscosity ratios in both 3-D and equivalent axisymmetric simulations. Tang et al. (Reference Tang, Adcock and Mostert2023) examined the influence of a range of
$\mathit{Oh}_d$
values less than
$0.1$
on the time of the start of breakup of the bag, discovering an exponential increase with
$\mathit{Oh}_d$
. All the aforementioned computational studies operate within
$\rho \approx 800$
, and thus, result in a weak dependence of the threshold to
$Oh_d$
, similar to experiments. This is apparent in figure 8 of Yang et al. (Reference Yang, Jia, Che, Sun and Wang2017), where the
$\mathit{We}_{cr}$
is observed to be within
$10\lt \mathit{We}_{cr}\lt 20$
for all
$\mathit{Oh}_d\lt 0.1$
. However, the threshold
$\mathit{We}_0$
can exhibit substantial variation with
$\mathit{Oh}_d$
for lower density ratios and ambient Ohnesorge numbers (
$\mathit{Oh}_o$
). This is highlighted by some (although limited) computational studies, that have explored the effect of
$\mathit{Oh}_d$
at lower density ratios. Han & Tryggvason (Reference Han and Tryggvason2001), for a density ratio of
$10$
, observed a significant decrease in the amount of deformation for higher drop viscosities. This decrease in deformation was speculated to lead to an increase in
$\mathit{We}_{cr}$
values. Kékesi et al. (Reference Kékesi, Amberg and Prahl Wittberg2014) linked the fragmentation morphology to the ratio of the characteristic times for shear and bag breakup, proportional to the ratio of viscosity ratio and density ratio, supported by 3-D simulations. Farsoiya et al. (Reference Farsoiya, Liu, Daiss, Fox and Deike2023) discovered that
$\mathit{We}_{cr}$
(based on turbulence dissipation rate) low density ratio (
$\rho =1$
) drops in isotropic turbulence are significantly influenced by the viscosity ratio.
Villermaux & Bossa (Reference Villermaux and Bossa2009) was the first to analytically describe the bag breakup process for an inviscid drop and derived a constant threshold value of 6 for
$\mathit{We}_{cr}$
, an underestimation compared with experimentally seen threshold values. Their work was extended to include the viscosity of the drop fluid, first by Kulkarni & Sojka (Reference Kulkarni and Sojka2014) and most recently by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021), resulting in a function of
$\mathit{Oh}_d$
that describes
$\mathit{We}_{cr}$
. This was corrected for the underestimation and led to a great match with previous experimental results. The analytical works mentioned above do not take into account ambient fluid properties such as ambient viscosity (
$\mathit{Oh}_o$
) and density (
$\rho$
) (which in turn dictates the relative velocity of the drop with the ambient) in influencing the resulting deformation characteristics. However, the density and viscosity contrasts between the ambient and drop fluids are generally very large for experimental systems, thus minimising the relative significance of ambient fluid properties relative to that of the drop. This allows for a good match between the analytical solution and corresponding experimental
$\mathit{We}_{cr}$
values, even with the aforementioned assumptions. However (as will be explored in detail in this work), for systems where the contrast between the (physical properties of) ambient and drop fluids is not substantial, these factors must be taken into account for the correct estimation of threshold
$\mathit{We}_0$
values.
The initial Reynolds number
${Re}_0$
(or alternatively the ambient Ohnesorge number
$\mathit{Oh}_o$
) also remains to be exhaustively explored, especially in the context of critical drop breakup threshold. Han & Tryggvason (Reference Han and Tryggvason2001) did simulations for different
$\mathit{Re}_0$
values for some low density-ratio cases, and observed a large reduction in drop deformations for low
$\mathit{Re}_0$
values. They speculated that this reduction in deformation might lead to an increase in
$\mathit{We}_{cr}$
values. Very few other works have explored or commented on the role of
$\mathit{Oh}_o$
on drop breakup (Guildenbecher et al. Reference Guildenbecher, López-Rivera and Sojka2009; Jain et al. Reference Jain, Tyagi, Prakash, Ravikrishna and Tomar2019; Marcotte & Zaleski Reference Marcotte and Zaleski2019). Jain et al. (Reference Jain, Tyagi, Prakash, Ravikrishna and Tomar2019) once again was one of the very few works to analyse the impact of
$\mathit{Re}_0$
on high density-ratio drops (
$\rho =1000$
) and observed higher incidences of plume formation in backward bag morphology for higher
$\mathit{Re}_0$
values.
Hence, there is the need for a single cohesive study analysing the effect of all the relevant non-dimensional parameters,
$\mathit{Oh}_o$
,
$\mathit{Oh}_d$
and
$\rho$
on drop deformation and breakup, and their impact on the threshold Weber number observed for critical breakup (
$\mathit{We}_{cr}$
). In this study we quantify the effect of each of these parameters using high-fidelity multiphase flow simulations. It should be noted that a distinction between liquid–gas and liquid–liquid drop-ambient systems has been maintained in the existing literature. However, fundamentally, the only differentiating factor between the two systems is the density and viscosity ratios, as well as the surface tension of the fluid interface. It is expected that the need for this distinction should vanish for a study that covers a sufficiently large parameter sweep involving
$\rho$
,
$\mathit{Oh}_o$
and
$\mathit{Oh}_d$
. Thus, a wide range of values of
$\rho$
and
$\mathit{Oh}_o$
have been quantified. To focus at a range rarely explored in existing studies,
$\mathit{Oh}_d$
values explored in this study are restricted to
$\mathit{Oh}_d \leqslant 0.1$
.
The paper starts with a description of the relevant impulsive acceleration problem (§ 2.1), the high-fidelity numerical model (§§ 2.2 and 2.3). The parameter space to be numerically explored is described in detail in § 2.4. The effects of
$\mathit{Oh}_d$
,
$\mathit{Oh}_o$
and
$\rho$
on the drop deformation, given that other parameters are constant, are described in detail, illustrating the forces and internal flow observed in the drops (§ 3). During the course of the parameter sweep, by simulating a range of Weber number values for every non-dimensional parameter set, the corresponding critical Weber number can be discovered. Based on the insights gained from the simulations, a novel non-dimensional parameter (
$C_{\textit{breakup}}$
) has been derived, incorporating the effects of all the relevant non-dimensional numbers. Here
$C_{\textit{breakup}}$
is found to be more effective in predicting the threshold of drop fragmentation, over the currently used
$\textit{We}_{cr}$
.
2. Problem description and formulation
2.1. Problem description and non-dimensionalisation
Let us consider a drop of diameter
$D$
containing a fluid of density
$\rho _d$
and dynamic viscosity
$\mu _d$
(subscript
$d$
implies properties associated with the drop). It is impulsively accelerated by a uniform flow of velocity
$V_0$
, density
$\rho _o$
and viscosity
$\mu _o$
(subscript
$o$
implies properties associated with the ambient medium, i.e. outside the drop), starting at
$t=0$
. The surface tension of the drop-ambient interface is
$\sigma$
.
Based on these initial conditions, we can define an initial Weber number
$\mathit{We}_0$
(
$\rho _o V_0^2 D/\sigma$
), which represents the competition between the dynamic pressure forces that drive the deformation of the drop and the capillary forces that resist this deformation. An increase in
$\mathit{We}_0$
corresponds to an increase in the maximum deformation achieved by the drop before it retracts to its equilibrium shape of a sphere, due to the action of surface tension. We then expect there to be a maximum
$\mathit{We}_0$
for which the surface tension forces barely prevent fragmentation in the drop. This threshold is called the critical Weber number
$\mathit{We}_{cr}$
(
$\rho _o V_{cr}^2 D/\sigma$
), where
$V_{cr}$
corresponds to the critical (lowest possible)
$V_0$
required to consistently realise a non-vibrational breakup. A Buckingham-Pi analysis ((2.1)) for this system reveals that
\begin{equation} \mathit{We}_{cr} = F \left ( \rho , \mathit{Oh}_o, \mathit{Oh}_d \right ). \end{equation}
The density ratio
$\rho$
= (
$\rho _d/\rho _o$
) is a measure of the inertia of the drop relative to the ambient medium, reflecting its responsiveness to external forces. The dynamic pressure forces exerted by the ambient medium scale with
$\rho _o$
, i.e. its effectiveness in inducing accelerations in specific parts (such as the peripheral rim or the core) or in the entirety of the drop is inversely proportional to
$\rho$
. The drop Ohnesorge number
$\mathit{Oh}_d =$
(
$\mu _d/\sqrt {\rho _d\,\sigma \,D}$
) is a ratio of capillary and momentum diffusion time scales and provides an estimate of how the energy supplied to a drop by external forcing is distributed among surface energy and viscous dissipation. The ambient Ohnesorge number
$\mathit{Oh}_o =$
(
$\mu _o/\sqrt {\rho _o\,\sigma \,D}$
) provides a non-dimensional, velocity-independent analogue for the initial Reynolds number
$\mathit{Re}_0$
(since
$\mathit{Oh}_o = \sqrt {\mathit{We}_0}/\mathit{Re}_0$
), and represents ambient viscosity given that other parameters are the same. Finally, we define an instantaneous Reynolds number
$\mathit{Re} =$
(
$\rho _0 V_{\mathit{rel}} D_{\mathit{rel}}/\mu _o$
) in addition to
$\mathit{Re}_0$
, in order to accommodate the generally significant change in centre-of-mass accelerations and frontal area across the time scales over which a drop deforms and fragments. Here
$\mathit{Re}_0$
is based on the drop’s instantaneous velocity deficit (
$V_{\mathit{rel}}$
) with the ambient medium and its frontal radius for the instantaneous deformed shape (
$D_{\mathit{rel}}$
).
If we choose
$\rho _o$
,
$V_0$
and
$D$
as the basis variables for non-dimensionalising the system, the dimensional variables can be non-dimensionalised as
$\tilde {\rho _d} = \rho$
,
$\tilde {\sigma } = 1/\mathit{We}_0$
,
$\tilde {\mu }_o = \mathit{Oh}_o/ \sqrt {\mathit{We}_0}$
and
$\tilde {\mu }_d = \mathit{Oh}_d \sqrt {\rho /\mathit{We}_0}$
.
Hence, for a specific ambient-drop fluid combination (i.e. fixed physical properties), a specific drop diameter and a specific inflow velocity, a set of
$\{ \rho , \mathit{Oh}_o, \mathit{Oh}_d, \mathit{We}_0 \}$
can be obtained that fully characterises the system. A drop can be simulated in Basilisk to assess whether it undergoes a non-vibrational breakup for the corresponding non-dimensional set. If the drop does not fragment,
$\mathit{We}_0$
is systematically increased (attributable to a decrease in
$\sigma$
in non-dimensional space or an increase in inflow velocity in dimensional space) and the impulsive acceleration simulation is rerun. These steps are repeated until the drop exhibits a non-vibrational breakup. The corresponding minimum
$\mathit{We}_0$
that marks the onset of non-vibrational breakup is the critical Weber number
$\mathit{We}_{cr}$
for that particular non-dimensional set
$\{ \rho , \mathit{Oh}_o, \mathit{Oh}_d, \mathit{We}_0 \}$
.
2.2. Model description
The simulations have been performed using the open-source solver Basilisk, which is well validated for two-phase flows on adaptive Cartesian meshes across a wide range of densities and viscosity ratios (Popinet Reference Popinet2003, Reference Popinet2009; Marcotte & Zaleski Reference Marcotte and Zaleski2019). For this work, we employ its two-phase Navier–Stokes solver. The numerical scheme is detailed in Appendix A; this section outlines only the assumptions and parameters specific to our problem.
The extensive parametric sweep required by this study (2.1) makes 3-D fragmentation simulations computationally infeasible. We therefore utilise axisymmetric simulations, since the primary deformation process – from initial flattening to the onset of bag inflation – is predominantly axisymmetric. We note that asymmetries develop at later stages (
$t^*\gt 1$
) due to more prominent instabilities in the interface and ambient flow, making a quantitative point-by-point accurate prediction using axisymmetric simulations difficult. However, the goal of this work is to characterise the drop’s deformation pathway and ultimate fate by analysing its qualitative temporal evolution, such as periods of growth/decay of the drop dimensions and the concavity of the corresponding temporal evolution. We hypothesise that such predictions do not require point-by-point quantitative accuracy. Therefore, axisymmetric simulations are deemed sufficient for this study. The validity of this assumption is assessed in § 2.3.
The axisymmetric domain used for all simulations in this work. At
$t=0$
, the simulation starts with a stationary axisymmetric spherical drop under impulsive acceleration.
The general simulation domain, used for the simulations in this study, is illustrated in figure 2. It is a square (for compatibility with quadtree meshes) axisymmetric domain of size
$L$
and its dimensions are carefully chosen to ensure that the drop always remains a sufficient distance from the boundaries. Here
$L$
can be as small as
$16$
for drops with high inertia (
$\rho \geqslant 100$
) and as large as
$64$
for drops with low inertia (
$\rho = 10$
). The top boundary is symmetric (
$\partial _nP=\partial _nu_t=u_n=0$
) and the bottom boundary is the axisymmetry axis. The left boundary allows a uniform ambient fluid inflow into the domain (
$u_n = V_0 = 1$
) and the right boundary allows the flow to exit the domain freely (
$\partial _nu_n=0$
). A drop of diameter
$D=1$
is initialised with its centre on the axisymmetric axis and its initial velocity set to
$0$
.
To ensure sufficient solution accuracy, we enforce the maximum allowed wavelet errors to
$\chi _u = 10^{-4}$
and
$\chi _c = 10^{-6}$
for velocity and volume fraction fields, respectively. The maximum allowed residual for the Poisson solve,
$\epsilon _p$
, is set to
$10^{-4}$
. The minimum allowed cell size is set to
$512$
cells per diameter, which corresponds to
$2^{13}$
elements (
$N=13$
levels of refinement) for
$L=16$
, except for cases with
$\mathit{Oh}_o=0.0001$
, which correspond to the highest values of
$\mathit{Re}_0$
, for which we use
$1024$
cells per diameter (
$N=14$
for
$L=16$
). Test simulations are run with different values for
$\chi _c$
,
$\chi _u$
,
$\epsilon _e$
and
$N$
, detailed in Appendix B, and the resulting volume fraction fields and streamwise and transverse lengths of the drops are compared. The results show that the chosen values of
$\chi _c$
,
$\chi _u$
,
$\epsilon _e$
and
$N$
are sufficient to ensure that the simulations result in drop shapes that are converged with respect to these parameters.
A significant fraction of the simulations in this study involve high density ratios (
$\rho \gt 500$
). For such high density ratios, a sharp interface can induce instabilities at the interface due to an unnatural spike in kinetic energy (Jain et al. Reference Jain, Prakash, Tomar and Ravikrishna2015). The phenomenon has also been observed in all of our large
$\rho$
simulations with low
$\mathit{Oh}_d$
(
$\mathit{Oh}_d \leqslant 0.001$
), albeit not presented here. In these cases, the upstream face shows unnaturally large surface instabilities, which lead to the removal of micro-droplets from the main drop. To overcome these issues, the interface is smeared by vertex averaging the volume fraction field. This approach helps to reduce density gradients across the interface, preventing its premature breakup. The numerical scheme with this smearing will be validated for a high
$\rho$
case in the next section.
At t = 0, the ambient fluid is quiescent and the drop is initialised with zero initial velocity. Given the incompressible nature of the flow and an infinite propagation speed of any information across the domain, the end of the first time step sees the entire domain attain a flow consistent with the left inflow boundary conditions. This involves the establishment of an incompressible flow around the drop, requiring the velocity field to be solenoidal. In real life, this process occurs in a finite amount of time, dependent on the velocity of the acoustic wave velocity. However, in a numerical system, this occurs in a single time step and leads to a jump in the drop centre-of-mass velocity, without gaining any corresponding deformation. The magnitude of this velocity jump in its centre-of-mass velocity
$\delta V_{cm}$
is inversely proportional to
$\rho$
(Marcotte & Zaleski Reference Marcotte and Zaleski2019). Hence, the effective relative velocity at the undeformed state of the drop reduces to
$V_{\mathit{eff}}=1-\delta V_{cm}$
. Accounting for this effective velocity becomes imperative when calculating the associated
$\mathit{We}_0$
of the system. To address this, simulations have been conducted for each pertinent
$\rho -\mathit{Oh}_o$
pair, resulting in the following observed velocity jumps – a substantial jump of approximately
$0.14$
for
$\rho =10$
, a jump of
$0.030$
for
$\rho =50$
, and a negligible velocity jump of
$0.015$
$(\sim \! 1.5\,\%)$
for
$\rho =100$
. For all simulations with
$\rho \gt 100$
, this first time-step jump is deemed insignificant. All non-dimensional parameters in the current work have been corrected to incorporate this jump.
2.3. Verification of axisymmetric drop simulations
To verify the capabilities of Basilisk in modelling high density-ratio drops subjected to impulsive acceleration, we reproduce the Bag breakup experiment as described in Flock et al. (Reference Flock, Guildenbecher, Chen, Sojka and Bauer2012). An ethyl alcohol drop is released from a certain height above an approximately uniform jet of air. The drop undergoes a nearly quiescent free fall for a height of
$188$
mm before entering a jet of air of mean velocity
$10$
m s−1 and a peak velocity of
$15$
m s−1. The drop, having acquired some vertical velocity during its fall, has a shape that is close to but not perfectly spherical when it enters the air jet. The drop then deforms as a result of aerodynamic forces exerted by the air jet and finally breaks up according to a bag breakup morphology. As the drop enters the jet, it initially experiences aerodynamic forces applied by the boundary layer of the flow, and then moves into the main flow with peak flow velocities.
A simplified axisymmetric version of this experiment is simulated in Basilisk with non-dimensional parameters derived from the dimensional parameters specified in Flock et al. (Reference Flock, Guildenbecher, Chen, Sojka and Bauer2012). The simulation parameters are
$V=1$
,
$\mathit{Oh}_o=2.3 \times 10^{-3}$
,
$\mathit{Oh}_d=5.9652 \times 10^{-3}$
, the choice of
$\mathit{We}_0$
depends on the choice of air-jet velocity between
$10$
and
$15$
m s−1,
$L=16$
and
$D=1$
. The simulation differs from the experiment in multiple, although minor ways. Firstly, the initial free fall of the drop is omitted since a gravity force perpendicular to the jet direction would render the system non-axisymmetric. Consequently, the slight deformation of the drop just before encountering the air jet is captured in the simulation. Secondly, unlike the experiments where the drop passes through a boundary layer of thickness approximately
$3$
mm before experiencing the peak
$15$
m s−1 jet velocity, the simulations instantaneously load the drop with the full velocity of the air jet. Considering that the provided
$\mathit{We}_0=13$
is based on the mean jet velocity, it will be essential to find the
$\mathit{We}_0$
appropriate for our simulation conditions (instantaneous loading), corresponding to velocities between
$10$
and
$15$
m s−1.
Panels (a) and (b) compare the streamwise and transverse lengths (see figure 4 for definitions of
$e_x$
and
$e_r$
) obtained from the simulation to various experiments: JA (2021) (Jackiw & Ashgriz Reference Jackiw and Ashgriz2021), Flock et al. (Reference Flock, Guildenbecher, Chen, Sojka and Bauer2012), and axisymmetric and 3-D simulation results from LM (2023) (Ling & Mahmood Reference Ling and Mahmood2023).
Figure 3 plots the streamwise (
$e_x$
) and transverse (
$e_r$
) lengths of the drop as a function of time, obtained through an axisymmetric simulation in Basilisk, and compares them to the experimental and simulation results from some recent works. The streamwise and transverse lengths obtained from the experimental data in Flock et al. (Reference Flock, Guildenbecher, Chen, Sojka and Bauer2012) is plotted, along with the experiments performed by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021). Lengths
$e_x$
and
$e_r$
corresponding to axisymmetric and 3-D simulations performed by Ling & Mahmood (Reference Ling and Mahmood2023) are also plotted for reference. The experiments performed by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021) had
$\mathit{Oh}_o$
and
$\mathit{Oh}_d$
values close to Flock et al. (Reference Flock, Guildenbecher, Chen, Sojka and Bauer2012), and hence, are a good reference for comparison.
It is observed that the current axisymmetric simulations reasonably capture the streamwise length until
$t^*\approx 1.2$
and transverse length until
$t^*\approx 1.5$
. The prediction of both the parameters worsens beyond this point. This coincides with the period of rapid streamwise and transverse expansion that follows the initial inflation of the drop into a bag-like structure. This is evident in the temporal evolution of the transverse length,
$e_r$
, as shown in figure 3(b), which exhibits a rapid growth in the drop size after
$t^*\approx 1.5$
. The constraint of azimuthal symmetry, inherent to axisymmetric simulations, limits the accurate representation of the intrinsically 3-D phenomena associated with bag inflation and subsequent rupture. These phenomena are driven by surface instabilities (Lozano et al. Reference Lozano, García-Olivares and Dopazo1998; Bremond & Villermaux Reference Bremond and Villermaux2005; Villermaux Reference Villermaux2007; Zhao et al. Reference Zhao, Liu, Cao, Li and Xu2011). In addition, axisymmetric simulations tend to overestimate the stagnation pressures at the downstream pole of the drop, which artificially restricts the inflation of the bag (Ling & Mahmood Reference Ling and Mahmood2023). In contrast, 3-D simulations do not suffer from such limitations. The 3-D simulations performed by Ling & Mahmood (Reference Ling and Mahmood2023) show a substantially better agreement with experiments beyond
$t^*\approx 1.2$
. However, the current axisymmetric simulations exhibit excellent agreement with the axisymmetric results reported by Ling & Mahmood (Reference Ling and Mahmood2023), which employed the same boundary conditions and numerical methods. Therefore, the numerical set-up employed in this work can be considered comparable to axisymmetric simulations in the recent literature.
It is worthwhile to note that the 3-D simulations show a slight overestimation in streamwise deformation compared with the experiments. This manifests as a lower minimum
$e_x$
in the 3-D simulations compared with the experiments. This
$e_x$
value corresponds to the drop’s streamwise dimension at the onset of bag expansion, called the initiation time by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021). This can be attributed to the difference in initial Weber number
$\mathit{We}_0$
between the experiments (
$\textit{We}_0 = 13$
) and the simulations (
$\textit{We}_0 = 15.3$
). This overestimation may also arise from the presence of gravity in the experimental set-up that leads to a non-spherical initial shape of the drop when it enters the air jet.
Even with the observed limitations of axisymmetric simulations in quantitatively predicting the point-by-point evolution of the drop dimensions, the simulations still capture the qualitative evolution of the drop morphology. The almost monotonic decrease in streamwise dimension until reaching an initiation time, and the subsequent growth resulting in a positive concavity, is captured by the axisymmetric simulations. In figure 3(a) the axisymmetric simulations agree reasonably well with properties such as initiation time, and the constant radial expansion rate after balancing time
$T_{bal}$
, i.e. the time when aerodynamic and surface tension forces balance (Jackiw & Ashgriz Reference Jackiw and Ashgriz2021). The axisymmetric simulations also succeed in achieving a similar initial pancake shape of the drop, as well as the subsequent appearance of a bag for such high density, low viscosity drops (Ling & Mahmood Reference Ling and Mahmood2023). For the purposes of this work, this is sufficient information to conduct a broad categorical analysis on the sensitivity of features to parameters such as
$\mathit{Oh}_o$
,
$\mathit{Oh}_d$
and
$\rho$
. Thus, the axisymmetric simulations are deemed sufficient for such water–air like systems.
However, in order to justify the use of axisymmetric simulations for the extensive parameter space relevant to this work, we would need to verify the ability of axisymmetric simulations in capturing the pancake shape and subsequent bag formation for a wider range of density ratios and ambient Ohnesorge numbers. This is done in Appendix C, where we compare 3-D and axisymmetric simulations for a few benchmark cases, with a special focus on non-trivial cases that show forward pancake and forward bag formation. In short, we observe an good agreement between 3-D and axisymmetric simulations for all cases except for
$\mathit{Oh}_o,\mathit{Oh}_d\leqslant 0.001$
, including the appearance of forward pancakes for high
$\mathit{Oh}_o$
cases and forward bags for low density-ratio cases. The case with
$\{\rho ,\mathit{Oh}_o,\mathit{Oh}_d\}=\{100,0.0001,0.001\}$
however shows a significant deviation in the post pancake deformation of the drop, with the axisymmetric simulations showing a forward bag formation while the 3-D simulations show a backward bag formation. Thus, any conclusions drawn from axisymmetric simulations for such low viscosity systems should be treated with extreme caution. For all other systems, we can be confident that the axisymmetric simulations capture the qualitative evolution of the drop morphology, including the appearance of forward pancakes and forward bags.
Hence, all information pertinent to this work can be reliably obtained through axisymmetric simulations. The axisymmetric simulations are significantly less expensive than 3-D simulations, and hence, well suited for the extensive parametric study proposed in this work.
2.4. Parameter space explored
The goal of this work is to systematically study the influence of the three non-dimensional parameters discussed, namely
$\rho$
,
$\mathit{Oh}_o$
and
$\mathit{Oh}_d$
, on drop deformation and breakup morphology. Table 1 summarises the parameter space explored through simulations in this work, encompassing 60 sets of
$\{\rho , \mathit{Oh}_o, \mathit{Oh}_d\}$
. Each set is simulated for different
$\mathit{We}_0$
values in order to identify its critical Weber number
$\mathit{We}_{cr}$
and the corresponding critical breakup morphology. This is achieved by simulating each
$\{\rho , \mathit{Oh}_o, \mathit{Oh}_d\}$
set with multiple
$\mathit{We}_0$
values to determine the lowest
$\mathit{We}_0$
value at which a non-vibrational breakup is observed, in short its
$\mathit{We}_{cr}$
.
A list of all values of
$\rho$
,
$\mathit{Oh}_o$
and
$\mathit{Oh}_d$
that form the part of the parametric space to be explored through simulations. In total, 60 sets of
$\{ \rho , \mathit{Oh}_o, \mathit{Oh}_d \}$
are considered, each is simulated for multiple
$\mathit{We}_0$
values to obtain
$\mathit{We}_{cr}$
. The minimum and maximum number of cells in the computational domain associated with all simulations is listed in the third column.
Given
$Re_0 \propto 1/\mathit{Oh}_o$
, a range of
$Re_0$
of the order of
$(10, 10^4)$
corresponds to a
$\mathit{Oh}_o$
range of
$(0.0001, 0.1)$
. Any higher
$\mathit{Re}_0$
becomes computationally infeasible due to significant turbulent vortices in the domain leading to a requirement for higher mesh resolution, lower time-step sizes and even 3-D simulations. Hence, this justifies the parameter space specified for
$\mathit{Oh}_o$
.
The influence of
$\mathit{Oh}_d$
on
$\mathit{We}_{cr}$
has been investigated extensively in the literature, primarily for systems with large
$\rho$
and
$\mathit{Oh}_o$
values, conditions typical of most experimental set-ups. For such systems,
$\mathit{We}_{cr}$
has been observed to be relatively constant (typically within
$10\lt \mathit{We}_{cr}\lt 20$
) for such systems for
$\mathit{Oh}_d\lt 0.1$
(Hsiang & Faeth Reference Hsiang and Faeth1995; Guildenbecher et al. Reference Guildenbecher, López-Rivera and Sojka2009; Yang et al. Reference Yang, Jia, Che, Sun and Wang2017). However, a similar comprehensive understanding of the influence of
$\mathit{Oh}_d$
is lacking for low
$\rho$
or
$\mathit{Oh}_o$
values. This study aims to address this gap by exploring the fragmentation threshold for the commonly seen values of
$\mathit{Oh}_d$
across a parameter space encompassing variations in both
$\rho$
and
$\mathit{Oh}_o$
beyond the typical experimental range. Naturally, to encompass the entire space of low and high density-ratio systems,
$\rho$
values are varied from
$10$
to
$1000$
.
The choice of such a comprehensive parameter exploration sets the stage for a detailed investigation into the nuanced interplay of these parameters on drop deformation and breakup characteristics.
3. Results
Non-dimensional pressure field
$P/(\rho _o V_0^2)$
renders for a drop undergoing backward bag fragmentation. Labelled A--E in the figure are the points of interest in a deforming drop and relate to the following features: the upstream pole/centre/core of the drop (A); the periphery of the drop (B); the rim of the drop, which in general has a higher local inertia compared with its centre (C); the downstream low pressure circulation zone, which can affect the motion of its rim if it is attached to the drop surface (unlike in the figure) (D); and the inflated bag, which inflates because of its low inertia and, hence, higher accelerations (E).
Figure 4 illustrates a typical drop deformation process, indicating the key features used to define several essential factors. Together, these factors provide an intuitive framework for understanding the deformation dynamics and will be used throughout this study.
-
(i) The local variation in inertia across the drop, which determines the local accelerations of the constituent parts, e.g. the difference in local inertia between the centre (A) and rim (C).
-
(ii) The pressure difference between the poles (A) and periphery (B), denoted by
${\Delta P}_{\textit{dri}v\textit{e}}$
.
${\Delta P}_{\textit{dri}v\textit{e}}$
is directly proportional to the stagnation pressures observed at the upstream pole of the drop. -
(iii) The surface stresses or viscous forces experienced by the upstream facing surface of the drop. These stresses are a function of the instantaneous Reynolds number
$\mathit{Re}$
, which can in many cases be approximated as equal to
$\mathit{Re}_0$
if the temporal evolution of centre-of-mass velocities and drop deformation do not significantly alter
$\mathit{Re}$
. -
(iv) The drop Ohnesorge number
$\mathit{Oh}_d$
, which dictates the distribution of the total energy supplied by the ambient flow to the drop between surface energy and the fluid momentum gained by the drop.
From the start of the deformation process until the formation of a distinct rim at
$t^* \approx 1$
, a drop does not exhibit any appreciable variation in local inertia along the lateral dimension. Hence, during this initial deformation phase, local inertia differences do not significantly influence the initial deformation; instead, the interplay between pressure and shear forces predominantly governs the process. Acceleration and, hence, the increase in velocity of the centre of mass is inversely proportional to total inertia and directly affects the instantaneous Reynolds number
$\mathit{Re}$
of the ambient flow past the drop. While
$\mathit{Re}$
dictates the shear stresses on the upstream surface, the relative velocity of the drop dictates the stagnation pressures at the upstream pole and, hence,
${\Delta P}_{\textit{dri}v\textit{e}}$
. The resulting pancake shape depends on the comparative strengths of the pressure difference and shear forces.
Once a drop develops local inertia variations across the lateral dimension as it deforms past the pancake stage, any subsequent deformation becomes strongly dependent on these variations in local accelerations. For the same external forces, regions of the drop with larger inertia experience much lower accelerations, and hence, lag behind their lower inertia counterparts.
The Reynolds number of the ambient flow past the drop dictates the strength, time scales, length scales and location of the downstream vortices (Forouzi Feshalami et al. Reference Forouzi Feshalami, He, Scarano, Gan and Morton2022). Thus, it is essential that we consider the interaction of these vortices with the rim for different Reynolds numbers to understand the final shape. The sensitivity of the rim to these flow characteristics is decided almost solely by inertia relative to the ambient fluid, i.e.
$\rho$
. A large density-ratio drop is expected to exhibit very little sensitivity to downstream vortices, and vice versa.
If we consider the specific drop case shown in figure 4, the ratio of the spatial extent of the drop along the axisymmetric axis to the spatial extent along the
$r$
axis provides its aspect ratio
$A_{xr}$
. This parameter will be used in future sections to quantify the deformation shown by the drops for the parameter space. In the first image (from left to right), we observe a flat pancake, which occurs when
${\Delta P}_{\textit{dri}v\textit{e}}$
predominantly drives the internal flow in the drop (over shear stresses). We also observe a clear toroidal rim (second image) that has a large local inertia and is therefore expected to lag the lower inertia centre of the drop. Due to the large inertia, the rim remains unaffected by the low-pressure zone created by the downstream vortex, which sheds a sufficient distance away from the rim and is not attached to the drop. Ultimately, the drop deforms into a backward bag morphology, where the centre inflates into a bag under the action of pressure forces at the stagnation point.
3.1. Density ratio
Panel (a) shows the temporal variation of (streamwise
$e_x$
and transverse
$e_r$
) axis lengths, and the
$x$
component of the centre-of-mass (cm) velocity of drops with different
$\rho$
values. The analytical relationship for
$\dot e_r$
as given by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021) is plotted for reference (see line labelled ‘JA(2021)’). Dashed lines represent linear fit lines of
$e_r$
from
$t^*=0.3$
to
$t^*=1.2$
. Internal velocity fields for
$\rho =10$
,
$\rho =100$
and
$\rho =1000$
are plotted in (b), (c) and (d), respectively. The upper half shows rvelocities (
$u_r$
), whereas the lower half of each plot shows
$x$
velocities (
$u_x$
). All drop systems presented have
$\mathit{Oh}_o=0.001$
,
$\mathit{Oh}_d=0.1$
,
$\mathit{We}_0=20$
.
This section illustrates the role of density ratio in influencing drop deformation and breakup morphology. Figure 5 illustrates the variation of the streamwise (
$e_x$
) and transverse (
$e_r$
) lengths of the drop and
$x$
-component of centre-of-mass (cm) velocities (
$V_{cm,x}$
), in the top and bottom of the plot, respectively. Density ratios from
$10$
to
$1000$
are shown. All cases show a decrease in their transverse length until
$e_x$
reaches a minima. We define this instant at which the drop achieves its minimum streamwise length as its pancake state. If the streamwise length of the pancake reaches a critical minimum, the drop unstably loses fluid from its core to its periphery, resulting in the formation of a thin fluid sheet near the core and a toroidal rim near the periphery, which eventually ruptures. The cases with a high density ratio, i.e.
$\rho = 100,500,1000$
, achieve their minimum
$e_x$
at
$t^* \approx 1$
, which coincides with the initiation time as described by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021), i.e. the time of initiation of inflation of the bag. Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021) suggested an initiation time of approximately
$1$
based on their experiments of a water drop in an air jet, which is consistent with the results presented here. All the high
$\rho$
cases also achieve a minimum
$e_x$
at
$t^* \approx 1$
. In addition, the minimum
$e_x$
, which can alternatively be interpreted as a maximum rim curvature, is nearly identical. Since all cases share the same surface tension coefficient (for the same Weber number), a likely hypothesis is that the drops deform streamwise until their surface tension forces, which are directly proportional to rim curvature, reach a similar magnitude. The rim curvature can be closely approximated by the reciprocal of the thickness of the flattened ‘pancake’ shape (Villermaux & Bossa Reference Villermaux and Bossa2009; Kulkarni & Sojka Reference Kulkarni and Sojka2014). Thus, for similar aerodynamic forces, the rim curvature is expected to be similar for two drops with the same surface tension coefficient.
In contrast to the high
$\rho$
cases, the low
$\rho$
cases (
$\rho =10,50$
) do not fragment for a initial Weber number of
$\mathit{We}_0 = 20$
. This aligns with the fact that the drops for the two cases do not achieve a minimum
$e_x$
as low as the high
$\rho$
cases. Instead of unstably losing fluid from the core to the periphery, these drops begin retracting back towards a spherical shape after reaching their maximum deformation or pancake state. This hints at the lower aerodynamic forces experienced by the low
$\rho$
cases, which dictates the magnitude of the corresponding surface tension forces required to balance the aerodynamic forces. This hypothesis is supported by examining
$V_{cm,x}$
of the lowest
$\rho$
drops, which rapidly lose their relative velocity with the ambient, and hence, experience lower aerodynamic forces at later times.
It is also observed that the temporal development of transverse length
$e_r$
very closely follows a linear trend after reaching the balancing time
$T_{bal}$
, which is consistent with the analytical relationship for
$\dot e_r$
as given by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021). The analytical relationship for
$\dot e_r$
is given by
\begin{equation} \dot e_r = \left ( \dfrac {a}{2} \right )^2 \left ( 1 - \dfrac {128}{a^2 \mathit{We}} \right ) T_{bal}, \end{equation}
where
$a=6$
and
$T_{bal} = 1/8$
as specified by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021). The analytical relationship for
$\dot e_r$
is plotted in figure 5(a) as a dash-dotted line labelled ‘JA(2021)’ for reference. The linear fit lines for
$e_r$
, from
$t^*=0.3$
to the specific initiation time for each case are plotted as dashed lines and show a good agreement with (3.1).
Pressure field plots for drops with three different density ratios are plotted: (a)
$\rho =10$
, (b)
$\rho =100$
, (c)
$\rho =1000$
for a low
$\mathit{Oh}_o$
system. All drops referred to here have the following common parameters:
$\mathit{Oh}_o=0.001$
,
$\mathit{Oh}_d=0.1$
,
$\mathit{We}_0=20$
.
Figures 5(b,c,d) and 6(a,b,c), compare the internal velocity and pressure fields for drops of different density ratios for
$\mathit{Oh}_o = 0.001$
(
$\mathit{Re}_0 \approx 4472$
), respectively. Such high
$\mathit{Re}_0$
results in low shear stresses acting on the drop’s upstream surface, and almost the entirety of aerodynamic forcing is exerted through pressure forces. Until the formation of a distinct toroidal rim at
$t\approx 1$
, local inertia variations across the lateral dimension of the drop remain small and, therefore, do not significantly influence the corresponding local accelerations experienced by each region of the drop. Hence, during the pancake formation stage, the deformation is primarily dependent on the balance between the shear and pressure forces applied by the ambient medium.
Low
$\rho$
drops, having lower total inertias, experience larger centroid accelerations. This leads to lower relative velocities, instantaneous Reynolds numbers (
$\mathit{Re}$
) and stagnation pressures compared with high
$\rho$
drops under identical conditions. Thus, low
$\rho$
drops show a higher local and volume-averaged responsiveness to external forces (such as attached vortical structures). Consequently, the deformation process becomes highly sensitive to
$\rho$
. For instance, the pressure field for three different
$\rho$
values is presented in figure 6, and shows the lowest driving pressure forces (
${\Delta P}_{\textit{dri}v\textit{e}}$
) for
$\rho =10$
case, allowing even the low shear stresses (due to high
$\mathit{Re}_0$
) acting at the upstream surface of the drop to be a significant contributor to the initial deformation and internal flow of the drop. As
$\rho$
increases, the stagnation pressure (and hence
${\Delta P}_{\textit{dri}v\textit{e}}$
) grows, making the contribution of shear stresses increasingly irrelevant. Thus, for the
$\rho =1000$
drop, the internal flow is expected to be completely driven by the pressure forces. The internal flow plots shown in figure 5 confirm this behaviour. The internal flow for the
$\rho =10$
drop is highest near the upstream surface and decreases to nearly zero at the downstream pole (see
$u_x$
), with the corresponding velocity gradient pointing in a direction normal to the upstream surface. The internal flow for
$\rho =1000$
(shown in panel d) on the contrary, has the highest
$x$
velocities at the upstream pole and not at the periphery. All cases also show negligible lateral inertia differences across the drop before
$t^*\lt 1$
. Thus, for the
$\rho =10$
case, the high shear at the periphery leads to a greater local acceleration there than at the centre. We also note that the the lower relative velocity of the low
$\rho$
case with the ambient leads to even lower
$\mathit{Re}$
values, resulting in a vortex that is not fully detached from its periphery, increasing the induced stresses. This creates a differential acceleration between the pole and the periphery for the
$\rho =10$
case, resulting in a forward-facing pancake. The
$\rho =1000$
case, on the other hand, forms a flat pancake. Thus, the competition between shear stresses and
${\Delta P}_{\textit{dri}v\textit{e}}$
controls the orientation of the pancake.
Beyond the pancake stage (at
$t^* \approx 1$
), a prominent toroidal rim forms in all cases (figure 6). This concentrates mass at the periphery, increasing the rim’s local inertia relative to the drop’s thinning core. The subsequent evolution is dictated by the interplay between this variation in local inertia, aerodynamic forces and the downstream vortex structure, which varies significantly with
$\rho$
. However, for all three cases, the rim eventually gains enough mass to decelerate relative to the centre of mass, forming a bag-like structure of backward orientation. The rate of evacuation of the core fluid towards the periphery and the lateral stretching of the pancake is significantly different for the three cases, leading to different bag morphologies.
The
$\rho =10$
drop experiences the smallest pressure forces driving its internal flow, and hence, its core experiences the smallest rates of evacuation. It hence takes longer for a prominent rim to form, coinciding with a delayed flipping of the thinned pancake from forward to backward. The downstream vortex shed from the periphery is also the weakest for this case and cannot induce any significant stretching of the rim. The end result is substantially less deformation and no fragmentation in this case. Conversely, the
$\rho =1000$
drop experiences the largest
${\Delta P}_{\textit{dri}v\textit{e}}$
driving its internal flow, and hence, the rate of evacuation of the core fluid is very high. The high inertia makes the rim fairly insensitive to downstream vortices and large
$\mathit{Re}$
leads to the vortex detaching from the periphery early and cleanly. It is only the intermediate
$\rho =100$
case that shows a sufficiently strong downstream vortex that sheds close to the rim and produces larger lateral stretching rates, reflected in the larger transverse internal velocities
$u_r$
compared with the other two cases (figure 5). In fact, only an intermediate
$\rho$
allows large enough
$\mathit{Re}$
to generate a strong downstream vortex and yet low enough inertia for the rim to be sensitive to such forces. This stronger stretching coupled with the intermediate rates of evacuation of the core results in a backward bag with a core that has not completely evacuated, leading to a bag-plume morphology. A similar explanation for the formation of a plume is provided by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021), where a faster bag inflation rate (due to higher
$\mathit{We}_0$
in the paper) compared with the movement of drop fluid from the centre to rim leads to the presence of an undeformed core at the centre of the drop. The volume of fluid contained in this undeformed core dictates the specific breakup observed – bag-plume, multi-bag or sheet thinning. This deformation process is also shown in Marcotte & Zaleski (Reference Marcotte and Zaleski2019, figure 4) for a low
$\mathit{Oh}_o$
case where a variation in
$\rho$
from
$10$
to
$2000$
is accompanied with a shift in pancake and breakup morphology exactly as observed here.
Pressure field plots are plotted for two different
$\rho$
values for a high
$\mathit{Oh}_o$
ambient flow. Results are shown for (a)
$\rho =10$
and (b)
$\rho =1000$
. All drops shown here have the following common parameters:
$\mathit{Oh}_o=0.1$
,
$\mathit{Oh}_d=0.1$
,
$\mathit{We}_0=20$
.
Let us now shift our attention to the effect of
$\rho$
on high
$\mathit{Oh}_o$
cases. Figure 7 shows the evolution of the pressure field around the drops with time. In both cases,
$\mathit{Oh}_o$
is
$0.1$
, which corresponds to
$\mathit{Re}_0=44.72$
. For such a low
$\mathit{Re}_0$
, we expect the shear stresses on the upstream surface to be substantial. Hence, despite observing substantially higher upstream stagnation pressures (and hence higher
${\Delta P}_{\textit{dri}v\textit{e}}$
) for (b) (
$\rho =1000$
) compared with (a) (
$\rho =10$
),
${\Delta P}_{\textit{dri}v\textit{e}}$
still does not dominate over the shear stresses during pancake formation. As expected, we see a forward pancake at
$t^*=0.948$
for both cases. Additionally, because of the low
$\mathit{Re}_0$
of the flow, the ambient flow remains attached to the drop’s surface, eliminating the formation of any downstream circulation zones. Subsequently, as the drop core is evacuated and a distinct rim is formed, a backward bag remains the only possible morphology. Panel (b) shows a backward bag breakup, while panel (a) shows much smaller deformations and does not fragment, the lower deformation is consistent with lower relative velocities, resulting in lower external forces.
In conclusion, the morphology of the pancake depends on the competition between the pressure difference between the poles and the periphery of the drop, with the shear stresses acting on the upstream surface. A flat pancake is observed when
${\Delta P}_{\textit{dri}v\textit{e}}$
is dominant, whereas a forward-facing pancake is observed when shear stresses are dominant. As the drop deforms past the pancake stage, it forms a bag, which can be forward or backward, depending on the local inertia of the rim and the strength of the downstream vortex. Local inertia depends on
$\rho$
and the rate of evacuation of fluid from the drop core. The strength and location of downstream vortices, on the other hand, depend on
$\mathit{Oh}_o$
and the instantaneous velocity of the drop, which again depends on inertia
$\rho$
. Finally, under conditions where the drop experiences a higher rate of lateral stretching (dependent on
$\mathit{Re}$
) compared with the rate of evacuation of the core (dependent on
${\Delta P}_{\textit{dri}v\textit{e}}$
), we may also observe a plume.
3.2. Ambient Ohnesorge number
Panel (a) shows the temporal evolution of the (streamwise
$e_x$
and transverse
$e_r$
) axis lengths for different
$\mathit{Oh}_o$
values. The analytical relationship for
$\dot e_r$
as given by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021) is plotted for reference (see line labelled ‘JA(2021)’). Dashed lines represent linear fit lines of
$e_r$
from
$t^*=0.3$
to
$t^*=1.2$
. Velocity fields are plotted for
$\mathit{Oh}_o=0.1$
and
$\mathit{Oh}_o=0.001$
in (b) and (c), respectively. For all plots,
$\rho =1000$
,
$\mathit{Oh}_d=0.1$
and
$\mathit{We}_0=20$
.
Figure 8(a) illustrates the temporal evolution of drop streamwise (
$e_x$
) and transverse (
$e_r$
) lengths and the
$x$
component of centre-of-mass velocities (
$V_{cm,x}$
), spanning
$\mathit{Oh}_o$
values from
$0.1$
(
$\mathit{Re}_0=44.72$
) to
$0.001$
(
$\mathit{Re}_0=4472$
) for a high density-ratio system. The discussions on the temporal evolution of
$e_x$
and
$e_r$
made in § 3.1 carry ideas that can be utilised to interpret the results shown in figure 8(a). The two lower
$\mathit{Oh}_o$
cases show a decrease in
$e_x$
to identical minimum values, achieved at the initiation time. This decrease in
$e_x$
is accompanied by a linear increase in
$e_r$
after
$T_{bal}$
, the rate of increase very similar to the analytical relationship (3.1) for
$\dot e_r$
as given by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021). However, in the case of
$\mathit{Oh}_o=0.1$
the streamwise length shows two distinct minima. The first minima occurs when the drop starts forming a forward pancake at
$t^* \approx 0.6$
, as seen in figures 8(b) and 9(a). At
$t^* \approx 1.2$
, the forward pancake starts flipping as the drop forms a backward bag, leading to a second minima in
$e_x$
right at the onset of bag inflation, which coincides with the initiation time as defined by Jackiw & Ashgriz (Reference Jackiw and Ashgriz2021). The second minima reached by the drop in the
$\mathit{Oh}_o=0.1$
case is higher than the minima reached by the other two cases, which hints at the lower rim curvatures achieved by the drop before bag inflation. It is also observed that this case follows the analytical linear increase in
$e_r$
only upto the first minima, after which the rate of increase in
$e_r$
decreases. This is expected as the analytical relationship was derived for experimental systems that typically form flat pancakes for high
$\mathit{Re}$
systems.
The drop with the highest value of
$\mathit{Oh}_o$
experiences the largest drag forces given the same relative velocity, the majority of which is imparted by the shear stresses acting on its surface, owing to the circular symmetry of the pressure field around a cylinder in low
$\mathit{Re}$
flows. This case thus shows the largest centre-of-mass acceleration, resulting in the highest
$V_{cm,x}$
throughout the deformation process. The lower relative velocities of this case result in lower stagnation pressures and, thus, lower
${\Delta P}_{\textit{dri}v\textit{e}}$
values. This hints at the cause for the lower rim curvatures achieved by the drop before bag inflation.
Pressure fields for drops with (a)
$\mathit{Oh}_o=0.1$
and (b)
$\mathit{Oh}_o=0.001$
. For both drops,
$\rho =1000$
,
$\mathit{Oh}_d=0.1$
and
$\mathit{We}_0=20$
.
Pressure fields for two different
$\mathit{Oh}_o$
values for
$\rho =1000$
have been plotted in figure 9. All non-dimensional parameters except
$\mathit{Oh}_o$
are the same for the two cases. For the drop in figure 9(a),
$\mathit{Oh}_o = 0.1$
, i.e.
$\mathit{Re}_0$
is very low, which corresponds to a large outside viscosity. Thus, the flow does not detach from the drop surface and leads to large viscous stresses on the upstream surface and, consequently, larger centre-of-mass velocities. The drop in figure 9(b), on the other hand, has a
$\mathit{Re}_0$
value that is 100 times larger, leading to much smaller shear stresses and, consequently, smaller centre-of-mass velocities. The larger velocities of the drop shown in panel (a) lead to smaller stagnation pressures and, consequently, smaller
${\Delta P}_{\textit{dri}v\textit{e}}$
compared with the case shown in panel (b). It should be noted that the instantaneous Reynolds number
$\mathit{Re}$
(instead of
$\mathit{Re}_0$
) would be a more accurate descriptor of the effective shear stresses experienced by a drop. However, this discrepancy between
$\mathit{Re}$
and
$\mathit{Re}_0$
(resulting from non-zero centre-of-mass velocities and frontal radius growths) plays a minor role in influencing shear stresses compared with the two orders of magnitude change in
$\mathit{Re}_0$
between the two cases. For the drop shown in figure 9(a), owing to lower stagnation pressures, the shear stresses acting on its upstream surface dictate its initial internal flow and resulting deformation. This is clearly demonstrated by the internal flow plot shown in figure 8(b) for the drop, where velocities are the highest at its upstream surface and decreases to zero at its downstream pole, coinciding with the location of the largest shear stresses applied by the ambient flow. This dominance of shear stresses results in the formation of a forward-facing pancake.
Conversely, for the higher
$\mathit{Re}$
system depicted in figure 9(b), shear stresses are significantly lower, which coupled with the larger
${\Delta P}_{\textit{dri}v\textit{e}}$
(compared with panel a) results in a pressure dominated internal flow. Consequently, the drop deforms into a flat pancake and the highest internal velocities occur at its upstream pole rather than at the periphery.
In summary, the orientation of the pancake, determined by the competition between
${\Delta P}_{\textit{dri}v\textit{e}}$
and shear stresses, is influenced by both
$\rho$
, due to its significant impact on
${\Delta P}_{\textit{dri}v\textit{e}}$
, and
$\mathit{Oh}_o$
, due to its significant impact on the
$\mathit{Re}$
, and consequently, the shear stresses exerted on the drop.
As the drops deform further, both cases develop a prominent rim, resulting in a large disparity in local inertia between the rim and the centre (prospective bag) of the drop. For the drop in figure 9(a), the extremely low
$\mathit{Re}$
results in an attached flow downstream of the drop, preventing the formation of a downstream vortex. In contrast, the drop in panel (b) with its large
$\mathit{Re}$
develops a downstream vortex, but the large local inertia (i.e.
$\rho$
) allows the vortex to detach early from the drop. Hence, for both cases, the rim does not experience any additional forces that can counteract the impact of the large local inertia (smaller local acceleration) of its rim. We observe that the drop in panel (a) flips orientation from a forward pancake to a backward bag, while the drop in panel (b) deforms from a flat pancake to a backward bag.
Panels (a) and (c) plot
$y$
velocities and the pressure field for a drop with
$\mathit{Oh}_o=0.01$
; whereas (b) and (d) plot
$y$
velocities and the pressure field for a drop with
$\mathit{Oh}_o=0.001$
. For both cases,
$\rho =50$
,
$\mathit{Oh}_d=0.1$
and
$\mathit{We}_0=20$
.
It is worthwhile to note that decreasing
$\mathit{Oh}_o$
motivates the formation of a plume. The reason for this can be discerned from the velocity field plots shown in figures 10(a) and 10(b). The lower
$\mathit{Oh}_o$
case develops a stronger downstream vortex due to its larger
$\mathit{Re}_0$
, as is evident from the larger
$r$
velocities at its rim. The common mediocre density ratio of
$\rho =50$
renders the toroidal rims of the drops more susceptible to lateral stretching due to their interaction with the downstream vortices. However, the induced drag on the lower
$\mathit{Oh}_o$
case is larger leading to comparatively higher rates of stretching. Ultimately, the larger
$\mathit{Oh}_o$
case fragments with a simple backward bag breakup morphology. The lower
$\mathit{Oh}_o$
case instead develops a plume that leads to the formation of an annular bag between the centre and the periphery, i.e. backward bag-plume morphology.
A case with much lower
$\mathit{Oh}_o$
(
${\lt } 0.001$
) would have a Reynolds number firmly in the free-shear regime (Forouzi Feshalami et al. Reference Forouzi Feshalami, He, Scarano, Gan and Morton2022), producing smaller, stronger and faster shedding vortices that form much closer to the drop periphery. Consequently, the drop in the higher
$\mathit{Re}_0$
case would experience even higher stretching rates from the poles to the periphery, as it is subjected to stronger induced drag forces due to better access to the low pressure zones downstream. This scenario when coupled with a low density-ratio drop could prevent the forward pancake from ever flipping, leading to a forward bag formation. This is indeed observed in the threshold fragmentation of the low
$\rho$
cases with
$\mathit{Oh}_o=\{0.001,0.0001\}$
and
$\mathit{Oh}_d=0.1$
, where the drop forms a forward bag without ever flipping its orientation. The pressure field for these cases is shown in Appendix C, and its validity justified through a comparison with analogous 3-D simulations.
In short, the rate of stretching, and consequently, the size of the resulting plume, increases with increasing proximity and strength of downstream vortices on the rim (due to an increase in
$\mathit{Re}_0$
). For low
$\mathit{We}_0$
simulations, forward bags are only observed for smaller inertia drops when
$\mathit{Oh}_o$
values are small, with
$\mathit{Re}$
preferably in the shear layer instability regime (
$1000\leqslant \mathit{Re}_0\leqslant 10^5$
).
3.3. Drop Ohnesorge number
We start our investigation into the role of the drop Ohnesorge number by studying cases with high density ratios and low ambient Ohnesorge numbers. In such drop-ambient systems, drops due to their large inertia maintain high relative velocities with the ambient flow, yielding large stagnation pressures (
${\Delta P}_{\textit{dri}v\textit{e}}$
) and large instantaneous Reynolds numbers (
$\mathit{Re}$
). The low shear stresses acting on its upstream surface (due to high
$\mathit{Re}$
) allow
${\Delta P}_{\textit{dri}v\textit{e}}$
to dictate its internal flow, resulting in the formation of a flat pancake. Given that all other non-dimensional parameters are the same across two cases, the case with the larger
$\mathit{Oh}_d$
has a larger drop viscosity
$\mu _d$
, and is expected to have lower internal flow velocities and circulations. Higher
$\mu _d$
also results in an exponential decrease in the incidence of surface instabilities (Fuster et al. Reference Fuster, Agbaglah, Josserand, Popinet and Zaleski2009). While an increase in surface tension decreases the wavelength of the fastest growing surface waves, an increase in
$\mu _d$
increases the length and time scales for which a capillary wave generated by an instability remains intact. As a drop accelerates and its relative velocity with the ambient flow decreases, the effective acceleration of the drop relative to the ambient medium decreases. Consequently, a surface wave that might have been unstable at the start of deformation, might be stable at more advanced stages of deformation, under the condition that the length scales and time scales of instabilities are large enough for the given
$\mu _d$
(Goodridge et al. Reference Goodridge, Shi, Hentschel and Lathrop1997).
Figure 11 illustrates the pressure fields for two such cases with
$\rho = 1000, \mathit{Oh}_o = 0.001$
and a varying
$\mathit{Oh}_d$
. Both cases show very similar (high) stagnation pressures and an
$\mathit{Re}_0$
of
$4472$
results in a well-defined downstream vortex detached from their peripheries, indicative of their large local and total inertia. As expected, at
$t^*\approx 1$
, both cases form a flat pancake and the onset of rim formation. However, for the drop in panel (b), we observe a high-pressure zone at its upstream pole that hints at the initiation of a plume. The internal flow field in figure 12(b) reveals an instability at the upstream pole of the drop, motivating a flow from its periphery to its upstream pole, hugging its upstream surface. The smaller viscosity of the drop in panel (b) (100 times smaller) facilitates the development of a prominent capillary instability at its pole and the corresponding pancake-plume shape. Furthermore, the location of the instability (i.e. the upstream pole) is a stagnation point and sees the highest accelerations out of all regions of the drop. This matches with the definition of Rayleigh–Taylor instabilities and might be the primary mechanism behind the development of a plume of this kind. According to Villermaux (Reference Villermaux2007), Jalaal & Mehravaran (Reference Jalaal and Mehravaran2014), an increase in density discontinuity motivates the formation of Rayleigh–Taylor instabilities. Our current simulations also show this behaviour, as only the cases with
$\rho = 500$
or
$1000$
and for the lowest
$\mathit{Oh}_d$
values form an unstable plume.
As the drops in the two cases continue to deform and more mass is transferred from their cores to their rims, substantial variations in local inertia (and hence local accelerations) start to develop. Notably, for the drop in figure 11(b), the plume has grown further and the drop now has two high local inertia regions – its core and its rim. The annular region connecting its core and rim has lower local inertia compared with both these regions and, hence, accelerates downstream relative to both, resulting in the growth of a bag between the plume and the rim. Ultimately, this annular bag fragments, culminating in a backward bag-plume morphology, as seen in panel (b).
Pressure fields around drops of different
$\mathit{Oh}_d$
values: (a)
$\mathit{Oh}_d=0.1$
and (b)
$\mathit{Oh}_d=0.001$
. All drops in these plots have
$\rho =1000$
,
$\mathit{Oh}_o =0.001$
,
$\mathit{We}_0=20$
.
Internal flows for two different
$\mathit{Oh}_d$
values: (a)
$\mathit{Oh}_d=0.1$
and (b)
$\mathit{Oh}_d=0.001$
. (c) Zoomed-in view of
$t^*=1.2649$
for (b). All drops in these plots have
$\rho =1000$
,
$\mathit{Oh}_o=0.001$
,
$\mathit{We}_0=20$
.
For the drop in figure 11(b), it should be noted that a reduction in
$\mathit{We}_0$
(while keeping other parameters constant) still results in the development of an instability-driven plume at its upstream pole, albeit of a smaller size. For instance, when
$\mathit{We}_0=16$
for the drop in figure 11(b), it does not deform enough to exhibit fragmentation (of the bag-plume kind). Consequently, solely decreasing
$\mathit{We}_0$
is not sufficient to shift the breakup morphology from a backward bag-plume to a simple bag breakup for the specific (
$\rho$
,
$\mathit{Oh}_d$
and
$\mathit{Oh}_o$
) set. This makes a backward bag-plume breakup the critical breakup morphology for this case – a feature of the physical properties of the system described by (
$\rho$
,
$\mathit{Oh}_d$
and
$\mathit{Oh}_o$
), and not just a function of boundary conditions (i.e. inflow velocity or
$\mathit{We}_0$
).
In contrast to the plume in figure 11(b) that originates from an instability due to lack of sufficient viscous damping, a decrease in
$\mathit{Oh}_d$
can also lead to a plume similar to that in figure 10(b). One such example is shown in figure 13.
The
$y$
velocities for two different
$\mathit{Oh}_d$
values: (a)
$\mathit{Oh}_d=0.1$
and (b)
$\mathit{Oh}_d=0.001$
. Both cases have
$\rho =100$
,
$\mathit{Oh}_o=0.001$
,
$\mathit{We}_0=13$
.
If we focus our attention on the specific case of drops,
$\mathit{Oh}_d$
can be interpreted as the ratio of the capillary time scale
$T_\sigma$
to the viscous time scale
$T_\mu$
, i.e.
$\mathit{Oh}_d = T_\sigma /T_\mu$
. The capillary time scale
$T_\sigma$
(
$= \sqrt {\rho _d D^3/\sigma }$
) is defined as the duration for a capillary wave of wavelength
$D$
to traverse a distance of
$D$
; while viscous time scale
$T_\mu$
(
$=\rho _d D^2/\mu _d$
) represents the duration for momentum to diffuse across the drop (Popinet Reference Popinet2009). A smaller
$\mathit{Oh}_d$
hence implies a relatively small
$T_\sigma$
compared with
$T_\mu$
, indicating that the information about interface deformation travels much faster than the rate of transfer of momentum to the drop fluid across the diameter. Hence, the downstream vortices could apply some induced drag on the drop rim, causing local acceleration relative to the core and subsequent deformation. However, this induced drag may not setup equivalent flow throughout the entire drop fluid. This is evident from the
$y$
-velocity plots shown in figure 13, where the lower
$\mathit{Oh}_d$
case (panel b) shows larger
$y$
velocities at the rim, indicating greater rates of stretching compared with the higher
$\mathit{Oh}_d$
case (a). Hence, the case in panel (b) due to its larger
$T_\mu$
results in a plume. It is essential to note that the drop in this case only shows a plume once it has started to form a bag, and the initial pancake at
$t^*\approx 1$
is flat. In contrast, the plume in figure 11(b) develops very early in the deformation process, right at the instant of formation of the pancake. Hence, the two types of plumes fundamentally differ in their formation mechanisms.
The fluid interface for three cases with
$\rho =1000$
,
$\mathit{Oh}_o =0.001$
,
$\mathit{We}_0=15$
, and (a)
$\mathit{Oh}_d=0.1$
, (b)
$\mathit{Oh}_d=0.01$
and (c)
$\mathit{Oh}_d=0.001$
.
According to much of the existing literature, if
$\mathit{Oh}_d \leqslant 0.1$
,
$Oh_d$
tends to have minimal impact on drop breakup mechanism. Therefore, for most studies, the choice of
$Oh_d$
is not focused upon, as long as it is ensured to be lower than
$0.1$
. However, the simulations conducted and analysed in this study do not corroborate with this understanding.
Another example emphasising the effect of
$\mathit{Oh}_d$
on drop deformation and breakup morphology is shown through drop interface plots in figure 14. In panel (a) the drop never achieves large enough deformation to undergo breakup. The drop in panel (b) on the other hand shows bag breakup for the same parameters except for
$Oh_d=0.01$
. The lower deformations achieved by the drop in panel (a) can be attributed to higher fluid viscosity, which provides resistance to internal flow and dissipates energy supplied by the ambient flow through surface forces. For
$\mathit{Oh}_d=0.001$
in panel (c), the breakup type shifts from a simple backward bag to a backward bag-plume breakup, with the plume formation driven by Rayleigh–Taylor instabilities. In short, a decrease in
$Oh_d$
is expected to reduce the required critical Weber number
$\mathit{We}_{cr}$
for a backward bag breakup. Hence, for the same
$We_0=15$
, we observe that the drop in panel (c) exhibits a backward bag-plume breakup, which is a multimode breakup that we expect to manifest at a
$\mathit{We}_0$
higher than that required for a simple backward bag.
4. Discussion
In this work, a parameter sweep using axisymmetric simulations was performed for multiple values of Weber number for every set of
$\{\rho , \mathit{Oh}_o, \mathit{Oh}_d \}$
possible in the parameter space defined in § 2.4. From this vast set of simulation data, we set out to achieve two primary objectives: (1) extract the influence of each involved non-dimensional parameter involved – specifically
$\rho$
,
$\mathit{Oh}_o$
and
$\mathit{Oh}_d$
– on drop pancake and breakup morphology; and (2) determine both Critical Weber number values as well as corresponding critical breakup morphologies for each unique combination of
$\{\rho , \mathit{Oh}_o, \mathit{Oh}_d \}$
in our parameter space. The first objective was addressed in detail in § 3. This section focuses on the second objective.
4.1. The threshold of impulsive drop breakup
Figure 15 shows the variation of critical Weber number (
$\mathit{We}_{cr}$
) against the drop Ohnesorge number (
$\mathit{Oh}_d$
) for drops of different density ratios (
$\rho$
) and the outside Ohnesorge numbers (
$\mathit{Oh}_o$
). Here
$\mathit{Oh}_d$
takes three different values in the parameter space:
$0.1$
,
$0.01$
and
$0.001$
. For every
$\mathit{Oh}_d$
, a
$\rho$
value is represented by a coloured vertical line that shows the range of
$\mathit{We}_{cr}$
values obtained due to variation in
$\mathit{Oh}_o$
. The lower
$\mathit{We}_{cr}$
values generally correspond to lower
$\mathit{Oh}_o$
values and vice versa. Therefore, for each
$\mathit{Oh}_d$
value in the plot, there exist five coloured vertical lines corresponding to the
$5$
$\rho$
values explored in the parametric sweep. It should be noted that each coloured line has been offset from its
$\mathit{Oh}_d$
value by a different amount for preventing overlaps with other
$\rho$
lines and, hence, improve clarity. All the cases are also explicitly marked with a uniquely shaped marker corresponding to each fragmentation morphology. Finally, all the experimental data for
$\mathit{We}_{cr}$
explored through this work is shown as a translucent area in the background of the plot, also shown as explicit markers in figure 1.
Plot of
$\mathit{We}_{cr}$
against
$\mathit{Oh}_d$
. Dependence on
$\mathit{Oh}_o$
is represented using vertical lines, with their vertical extent representing corresponding variation in
$\mathit{We}_{cr}$
. Dependence on
$\rho$
is shown through different coloured vertical lines (offset from its true
$x$
location to prevent overlaps with other lines) representing each
$\rho$
value in the parameter space. Specific markers are used to represent all critical breakup morphologies observed in the simulations. The experimental data for
$\mathit{We}_{cr}$
from figure 1 is shown as a translucent area in the background of the plot.
On the basis of all the simulation results and figure 15, the following conclusions emerge.
-
(i) Consistency with experimental data: all simulations with high
$\rho$
(
$\geqslant 500$
) values critically fragment at
$\mathit{We}_{cr}$
values that very closely match historical experimental works. This is expected given that most experimental studies have historically focused on impulsive fragmentation for water--air analogous systems, which have
$\rho \geqslant 500$
. Only cases with
$\rho \lt 500$
show any appreciable deviation from the experiments. -
(ii) Sensitivity to
$\mathit{Oh}_o$
: the lowest
$\rho$
cases exhibit the largest variations in critical Weber numbers with respect to changes in
$\mathit{Oh}_o$
for the high
$\mathit{Oh}_d$
cases, as indicated by the length of the vertical lines in figure 15. A drop with a large
$\rho$
(
${\geqslant}500$
) experiences larger relative velocities, making
${\Delta P}_{\textit{dri}v\textit{e}}$
the dominant factor driving its deformation in most cases. Even in cases with large external shear stresses on the drop (the largest
$\mathit{Oh}_o$
cases), once a clear rim has formed in the deformed drop, local inertia variations take over the deformation process. Only the lower inertia (
$\rho$
) drops show any appreciable sensitivity to changes in
$\mathit{Oh}_o$
and the corresponding differences in downstream vortical structures. An exception exists when
$\mathit{Oh}_o$
is in the free-shear regime (
$\mathit{Re}_0\gt 10\,000$
(§ 3.2), and even the drops with a high
$\rho$
respond strongly to the resulting downstream vortices, fragmenting with a forward bag morphology. -
(iii) Variations in fragmentation morphologies: the critical breakup morphology transitions from a forward bag (
$\rho =10$
) to a backward-transition (
$\rho =50,100$
) (see figure 16) to a backward and backward-plume bag with increase in
$\rho \geqslant 500$
. As
$\rho$
increases, the drop’s rim is expected to start lagging behind the drop core at some point during its deformation, when local inertia of its rim becomes substantially larger than that of its core. The shift in morphology from a forward to a backward bag is observed only for low
$\mathit{Oh}_o$
values, which produce downstream vortices that are strong enough to compensate for the larger local inertia of the rim. Alternatively, large
$\mathit{Oh}_o$
cases always exhibit backward bag breakup at critical conditions. This is due to the downstream vortices being weak or non-existent for low
$\mathit{Re}_0$
flows. -
(iv) Sensitivity to
$\mathit{Oh}_d$
: in addition to the drop Ohnesorge number’s role in controlling the sensitivity of
$\mathit{We}_{cr}$
to
$\mathit{Oh}_o$
, a decrease in
$\mathit{Oh}_d$
also affects the critical fragmentation morphology by motivating the formation of a plume. This is seen for the lowest
$\mathit{Oh}_d$
(
$0.01, 0.001$
) and
$\mathit{Oh}_o$
(
$0.001, 0.0001$
) values, and for the largest
$\rho$
(
$500, 1000$
) drops (star shaped marker in figure 15). Such drops show an unstable plume at the upstream poles of their flat pancakes, which are also locations of maximum acceleration in the drops and can be attributed to Rayleigh–Taylor instabilities. Both low viscosity and larger density ratios motivate the development of such instabilities (Villermaux Reference Villermaux2007; Guildenbecher et al. Reference Guildenbecher, López-Rivera and Sojka2009; Jalaal & Mehravaran Reference Jalaal and Mehravaran2014). Since the pancake for this non-dimensional set starts with a plume even for lower non-critical
$\mathit{We}_0$
values, a plain backward bag breakup can never manifest for such systems, and hence, a backward bag-plume breakup becomes its critical breakup morphology. The drop Ohnesorge number
$\mathit{Oh}_d$
also influences the rate of evacuation of the drop core, which ultimately results in a forward bag for all high
$\mathit{Oh}_o$
cases, when shear stresses drive pancake formation.
Panel (a) is a path diagram of all deformation paths a spherical drop under impulsive acceleration can take when breaking up critically. A spherical drop can deform into three types of pancakes, each of which can further deform into one of four breakup morphologies, the corresponding
$\rho ,\mathit{Oh}_o,\mathit{Oh}_d,\mathit{We}_0$
parameter space is shown in the phase diagram (b). The blue-highlighted region in (b) indicates cases with the lowest
$\mathit{Oh}_o$
and
$\mathit{Oh}_d$
values, for which the axisymmetric simulations may not be representative of reality, as discussed in Appendix C.
The conclusions drawn from figure 15 show that the accepted idea of critical Weber number being almost independent of drop Ohnesorge number for
$\mathit{Oh}_d\lt 0.1$
might not always hold, especially for systems that stray too far from properties analogous to the water-air system. Furthermore, the critical breakup morphology need not necessarily be a backward bag breakup. Backward bag-plume and forward bag morphologies can be the critical morphologies for certain low
$\rho$
and low
$\mathit{Oh}_o$
cases.
A path diagram, as illustrated in figure 16(a), can be drawn that summarises all the deformation paths a spherical drop might follow as it deforms to critical fragmentation under impulsive acceleration. A companion phase diagram (figure 16
b) provides the non-dimensional parameter space that results in one of the four observed fragmentation morphologies. It is worth highlighting that all the breakup paths provided in the diagram are for their respective critical conditions, and consequently, the
$\mathit{We}_{cr}$
values corresponding to different paths need not be the same.
For small
$\rho$
or high
$\mathit{Oh}_o$
values or both, shear stresses drive the internal flow, resulting in a forward pancake. The fate of this forward-facing orientation is then contingent on the balance between local inertia differences and the strength and proximity of downstream vortices to the rim. For systems with large
$\mathit{Oh}_o$
, irrespective of density ratio, downstream vortices either do not form or shed further downstream. The drop is not subjected to large lateral stretching rates, allowing the formation of a prominent toroidal rim. The resulting lateral inertia differences lead to the flipping of the bag from forward to backward orientation, eventually fragmenting with a backward bag morphology. On the other hand, for small
$\mathit{Oh}_o$
(and, hence, small
$\rho$
for a forward pancake) coupled with large
$\mathit{Oh}_d$
, the strong, fast and proximal downstream vortices generate large lateral stretching rates, preventing the formation of a prominent rim. As a result, the drop continues to hold its forward-facing orientation and breaks up with a forward bag morphology.
When
$\rho$
is large and
$\mathit{Oh}_o$
is small, the drop under critical conditions either deforms into a flat pancake or a flat pancake with a plume depending on whether
$\rho$
and
$\mathit{Oh}_o$
are at the extreme ends of the parameter space. The largest values of
$\rho$
and the smallest values of
$\mathit{Oh}_o$
and
$\mathit{Oh}_d$
lead to the appearance of a plume at the upstream pole of the flat pancake. It can be hypothesised that the low viscosities of both outside and drop fluids do not provide sufficient viscous dissipation to stabilise the jet ejected at the upstream pole (due to Rayleigh–Taylor instability) of the drop. From this pancake shape, the only possible critical breakup morphology is a backward bag-plume breakup. All other intermediate cases form a flat pancake that can form either a backward-transition breakup (for intermediate
$\rho$
values and low
$\mathit{Oh}_o$
values) or a backward bag (for all the remaining cases). A backward-transition breakup is a forward bag with a flipped rim, i.e. a drop that at its final moments gains enough inertia in its rim to start the bag flipping process. As expected, this is observed for intermediate
$\rho$
values where neither local inertia differences nor downstream vortices outright dominate the dynamics of the drop’s periphery.
The path diagram, together with the phase diagram, only informs of the types of pancakes and corresponding general breakup morphologies observed for a specific set of
$\{\rho , \mathit{Oh}_o, \mathit{Oh}_d\}$
under threshold conditions. For information on the lowest
$\mathit{We}$
required to achieve the corresponding non-vibrational breakup (i.e.
$\mathit{We}_{cr}$
), one may refer to figure 15.
4.2. Bag inflation characteristics
An understanding of the time scales involved for the inflation of bags during bag breakups is essential in the correct estimation of bag burst time scales and effective centre-of-mass velocities for bag breakups. Villermaux & Bossa (Reference Villermaux and Bossa2009) were the first to give an analytical description of bag inflation rates for backward bag morphology. They found that the amplitude of bag inflation increased exponentially with time with an exponent factor of two for an inviscid drop. Kulkarni & Sojka (Reference Kulkarni and Sojka2014) extended the work of Villermaux & Bossa (Reference Villermaux and Bossa2009) to include drop viscosity and numerically obtained a similar exponential relationship for bag inflation, now also dependent on
$\mathit{Oh}_d$
. However,
$\mathit{Oh}_d$
was found to have a very small impact on bag inflation rates. On the other hand,
$\mathit{We}_0$
has a more dramatic effect on the exponent factors governing inflation, with higher
$\mathit{We}_0$
showing faster inflation rates. The large variety of simulations in the current work across a large parameter space present an opportunity to explore the role of the relevant non-dimensional parameters on bag inflation. Before we proceed, it is essential to emphasise that the axisymmetric simulations conducted in this work cannot capture the bag inflation too far beyond the initiation time (Ling & Mahmood Reference Ling and Mahmood2023; Tang et al. Reference Tang, Adcock and Mostert2023). Hence, the bag inflation rates presented here are not useful in their absolute values, and only representative of the early stages of bag inflation. However, the differences in bag inflation rates for different cases, solely in the context of other axisymmetric simulations, still provides useful insights into the functional dependence of early bag inflation rates on these parameters.
Bag inflation
$\alpha (t^*)$
with time
$t^*$
is shown for some unique backward bag breakup cases. The non-dimensional parameter set for each case is of the form
$\{ \rho , \mathit{Oh}_o, \mathit{Oh}_d, \mathit{We}_0 \}$
. The cases plotted here include simple backward bags (
$\{500, 0.001,0.1,17\}$
,
$\{1000, 0.1,0.1,20\}$
$\{1000, 0.001,0.1,20\}$
) and backward-plume bags (
$\{100, 0.001,0.1,20\}$
,
$\{500, 0.001,0.1,20\}$
,
$\{1000, 0.001,0.001,20\}$
). In addition,
$\{1000, 0.1,0.1,20\}$
initially forms a forward pancake shape that then flips to a backward bag. Through this plot, the effect of
$\mathit{We}_0$
,
$\mathit{Oh}_o$
and
$\rho$
on bag inflation rates is highlighted.
With this caveat in mind, we highlight some key observations related to the bag inflation rates for a few simulations where backward bag breakup is observed. The evolution of bag inflation
$\alpha (t^*)$
with time is illustrated in figure 17, where
$\alpha (t^*)$
is equal to the horizontal extent of the drop
$e_x$
(as shown in figure 4). For every plot, the drop shows a decrease in
$\alpha (t^*)$
up until it reaches the end of the pancake stage at
$t^* \approx 1$
, beyond which the drop begins to inflate. The only exception is the
$\{1000, 0.1,0.1,20\}$
case, which initially shows an increase in
$\alpha (t^*)$
due to deforming to a forward pancake. The bag inflation stage only starts at
$t^* \approx 2.5$
for this case. Once firmly in the bag inflation stage, all cases plotted here show an exponential growth
$\alpha (t^*)$
with time with a specific exponential growth factor, marked in figure 17 as numbers alongside each plot.
From the plots, a number of key observations can be made. The
$\{1000, 0.001, 0.001, 20\}$
case, which has properties closest to a water--air drop-ambient system, shows an inflation growth factor of 1.96, which is very close to the value that was analytically found by Villermaux & Bossa (Reference Villermaux and Bossa2009), and was matched against a bag breakup experiment of a water drop. Even though this case is a backward-plume breakup, its inflation growth rate matches that of a simple bag and is not affected by the presence of the plume.
For the
$\{500, 0.001, 0.1\}$
case, out of the two Weber numbers shown in the plot
$(\mathit{We}_0 = 17,\, 20)$
, the higher
$\mathit{We}_0$
drop shows a faster bag inflation rate, indicated by the larger exponential growth of its
$\alpha (t^*)$
. This is an expected observation, since the overall deformation rate is also higher for a higher
$\mathit{We}_0$
case, owing to the lower surface tension forces relative to the dynamic pressure forces driving the drop’s deformation. Hence, the same amount of energy (supplied to the drop by the external forces) should lead to a larger change in the interface area as a result of a lower surface tension.
We also observe that the
$\rho =100$
case shows a dramatically lower inflation rate compared with analogous
$\rho =500$
or
$\rho =1000$
cases. This may be attributed to the lower relative velocities with the ambient flow observed for the low
$\rho$
case, which results in a decrease in its effective Weber number, and hence, lower bag inflation rates.
Drop and ambient viscosities appear to have an inverse effect on bag inflation rates. Among the three
$\rho =1000$
cases shown in the plot, the higher drop and ambient Ohnesorge number cases show the larger exponential growth rates of
$\alpha (t^*)$
.
More generally, figure 17 shows that
$t \approx \tau$
(or
$t^* \approx 1$
) is a good representative time scale for the start of the bag inflation process. This observation is consistent with the aspect ratio plots shown in figures 5 and 8, where the pancake formation stage almost always ended at
$t^*\approx 1$
.
If the aerodynamic forces acting on a drop are large enough to initiate bag inflation in its pancake, it is expected that the drop will go on to or be very close to fragmentation. Therefore, we hypothesise that the balance of all forces acting on a drop at this stage (aerodynamic forces driving its deformation, and capillary and viscous forces resisting it), i.e. at the end of its pancake formation and at the start of its bag inflation, is representative of the overall stability of the drop fragmentation process. We use this understanding to obtain a better estimate of the effective aerodynamic forces on the drop in the next section.
4.3. A parameter for prediction of breakup threshold
Most previous studies have characterised the threshold for impulsive breakup of spherical drops using a Weber number based on the initial relative velocity with ambient (
$\mathit{We}_0=\rho _o V_0^2D/\sigma$
). For all cases with properties analogous to the water--air system, i.e.
$\rho \gt 500$
,
$\mathit{Oh}_o\lt 0.01$
and
$\mathit{Oh}_d\lt 0.1$
, critical breakup occurs consistently at a critical Weber number of
$\mathit{We}_{cr} \approx 14$
. However, as has been discovered and exhaustively described through the simulation results in the previous sections (summarised succinctly in figures 15 and 16), different cases that stray away from the non-dimensional space described by water--air systems do not show the same threshold Weber number value. Substantially higher
$\mathit{We}_{cr}$
values are observed for cases with low density ratios and high Ohnesorge numbers. The Weber number
$\mathit{We}_0$
represents the ratio of the pressure forces applied on a drop surface by the ambient medium, based on its initial relative velocity to the surface tension acting against any change in surface energy. However, drop deformation also depends on the viscous forces applied by the surrounding flow, the inertia and, hence, the acceleration, and the viscous dissipation against internal fluid flow. These effects have been explained in detail in § 3 through simulations for varying
$\mathit{Oh}_o$
,
$\mathit{Oh}_d$
and
$\rho$
values, respectively. Hence,
$\mathit{We}_0$
does not capture the role of all the factors relevant in the drop deformation process. We aim to construct a new non-dimensional group that aggregates the effect of all the parameters, namely
$\mathit{We}_0$
,
$\mathit{Oh}_o$
,
$\mathit{Oh}_d$
and
$\rho$
, and shows a consistent critical value demarcating the threshold of breakup under impulsive acceleration for the complete parameter space explored in this study.
Let us assume that this new non-dimensional number, denoted by
$C_{\mathit{breakup}}$
, is a function of all the dimensional variables involved in the breakup process (4.1). There are three independent dimensions in the problem. Hence, through Buckingham-Pi analysis, we can obtain at most four independent non-dimensional numbers. The four non-dimensional numbers relevant to this problem have already been defined in § 2.1, namely
$\rho$
,
$\mathit{Oh}_o$
,
$\mathit{Oh}_d$
and
$\mathit{We}_0$
. It is expected for
$C_{\mathit{breakup}}$
to be described as
\begin{equation} C_{\mathit{breakup}} = f (V_0, D, \rho _o, \rho _d, \mu _o,\mu _d, \sigma ), \end{equation}
which when non-dimensionalised gives us
\begin{equation} C_{\mathit{breakup}} = f (\mathit{We}_0, \rho , \mathit{Oh}_o, \mathit{Oh}_d). \end{equation}
We derive
$C_{\mathit{breakup}}$
using a method analogous to that of Blackwell et al. (Reference Blackwell, Deetjen, Gaudio and Ewoldt2015, & Section IV). In their work, Blackwell et al., developed a non-dimensional group to describe the stick-splash behaviour of a drop impacting a surface. This group was based on the ratio of forces promoting splashing (inertial forces) to those promoting sticking (dissipative forces). Similarly, we define
$C_{\mathit{breakup}}$
based on the competition between forces driving fluid radially outward from the drop’s core and forces resisting this outward flow. We hypothesise that this competition determines the threshold for secondary fragmentation (the breakup of the primary drop into smaller droplets), where
$C_{\mathit{breakup}}$
is defined as
\begin{equation} C_{\mathit{breakup}} = \frac {\text{Forces driving drop deformation}}{\text{Forces resisting drop deformation}}. \end{equation}
An alternate version of figure 15 where instead of
$\mathit{We}_{cr}$
, the variation of
$C_{\mathit{breakup}}$
with respect to
$\rho$
,
$\mathit{Oh}_o$
and
$\mathit{Oh}_d$
is plotted. In addition to simulation data, all available experimental data in the relevant non-dimensional space for critical bag breakup (plotted in figure 1) is also plotted here for reference.
The forces will be evaluated at
$t^*\approx 1$
, which represents the end of the pancake stage and the start of the bag inflation process. The choice of
$t^*\approx 1$
as the critical time controlling the bag inflation process has been justified in § 4.2.
Let us assume that a spherical drop of diameter
$D_0$
at
$t^*=0$
deforms in to a disk of diameter
$D$
, bounded radially by a semi-circular ring of diameter
$W$
, and has an aspect ratio
$A_{xr}$
at
$t^*=1$
(figure 18
a). Mass conservation when applied to the drop gives us
$D \approx (2/A_{xr})^{(1/3)} D_0$
and
$W \approx (2/3) (D_0^3/D^2)$
, if we ignore the mass that the curved periphery holds.
All regions of the pancake within the flat disk region (
$r\lt D$
) have a very large local radius of curvature, and hence, contribute negligibly to surface tension forces locally. The periphery is the only region with any appreciable curvature equal to the sum of local longitudinal and azimuthal curvatures, i.e.
$(2/W + 2/D) = 2(1+A_{xr})/W$
. Therefore, crossing the drop interface at its periphery produces a pressure jump of
$2\sigma (1+A_{xr})/W$
according to the Young–Laplace equation that acts on the narrow cylinder between the semi-circular peripheral ring and the flat disk of area
$\pi D W$
at a radial distance
$D/2$
(point B in figure 18
a). The total contribution of surface tension can be estimated by calculating the force applied by the excess pressure on this area, given as
\begin{equation} F_\sigma = 2 \pi \sigma D\, (1+A_{xr}). \end{equation}
The surface tension force is directed against the movement of the drop fluid from the core to the periphery.
If the flow around the pancake is assumed to be a potential flow with a stagnation point at the upstream pole, the pressure at the upstream surface of the pancake at radial distance
$r$
from the longitudinal (axisymmetric) axis can be estimated (Villermaux & Bossa Reference Villermaux and Bossa2009; Jackiw & Ashgriz Reference Jackiw and Ashgriz2021) as
\begin{equation} p_o(r) = p_o(0) - \dfrac {\rho _o a^2 V_{\mathit{rel}}^2}{8 D_0^2} r^2. \end{equation}
Here,
$p_o(0) = 0.5 \rho _o V_{\mathit{rel}}^2$
is the stagnation pressure at the upstream pole. Thus, pressure is the maximum at the upstream pole and drops as we move radially away from the upstream stagnation point towards the periphery. The potential flow around a rigid body has a stretching factor
$a$
of
$6$
for a sphere and
$4/\pi$
for a flat disk, the latter being applicable in this work for the flow around a flat pancake. The true pressure inside the drop is a superposition of aerodynamic pressures and excess pressure due to surface tension. Since the local curvature of the drop interface is negligible for
$0\lt r\lt (D-0.5W)$
, the pressure inside the drop
$p_d(r)$
near its upstream surface in the disk region is almost equal to
$p_o(r)$
. We can thus calculate the independent contribution of the aerodynamic pressure forces in driving the evacuation of the drop core – by integrating the infinitesimal force due to
$p_o(r)$
acting on a cylindrical surface of radius
$r$
and width
$W$
across the radius of the pancake:
\begin{align} F_p &= \int _0^{0.5D} \left [ p_o(r)\,(2\pi r W) - p_o(r+{\rm d}r)\,(2\pi (r+{\rm d}r)W) \right ] \end{align}
\begin{align} &= - 2\pi W \int _0^{0.5D} \left [ p_o(r) + \dfrac {\partial p_o(r)}{\partial r}\,r \right ] {\rm d}r, \end{align}
\begin{align} F_p &= \dfrac {a^2}{48} \pi \rho _o V_{\mathit{rel}}^2 D_0\, D - \dfrac {1}{3} \pi \rho _o V_{\mathit{rel}}^2 \dfrac {D_0^3}{D} = F_{p,d} - F_{p,s}. \end{align}
Here
$F_{p,d}$
is the contribution of the radial pressure drop in the ambient flow in support of evacuation of the core, whereas
$F_{ps}$
is the contribution of stagnation pressure against it.
In addition to surface tension and external dynamic pressures, we must estimate the contribution of drop fluid viscosity in dissipating (part of) the kinetic energy of the internal flow. An estimate for the differential viscous dissipation power can be obtained by using the relation (Batchelor Reference Batchelor1967; Deville, Fischer & Mund Reference Deville, Fischer and Mund2002; Rimbert et al. Reference Rimbert, Castrillon Escobar, Meignen, Hadj-Achour and Gradeck2020)
\begin{equation} \varPhi = 2\mu _d\, \boldsymbol{d}:\boldsymbol{d} = 2\mu _d\; {\rm d}_{ij}{\rm d}_{ij}, \end{equation}
where
$\boldsymbol{d} = 0.5 [ \boldsymbol{\nabla }\boldsymbol{u} + (\boldsymbol{\nabla }\boldsymbol{u})^T ]$
is the rate of deformation tensor,
$\boldsymbol{d}:\boldsymbol{d}$
is its dyadic product and
$\boldsymbol{u}$
is the velocity vector field. Using mass conservation, Kulkarni & Sojka (Reference Kulkarni and Sojka2014) derived a relation
$u_r = r\dot {D}/D$
for the radial velocity of the drop fluid in terms of its radial expansion rate, which in Cartesian coordinates can be written as
$\boldsymbol{u} = (\dot {D}/D)(x \boldsymbol{\hat {j}} + y\boldsymbol{\hat {k}})$
, given that the transverse plane is synonymous with the
$yz$
plane. Substituting this velocity function in (4.8) and dividing by the local velocity scale, we can obtain a scale for the small local viscous force, which is always oriented in a direction opposite to the internal flow velocity (i.e. negative radial direction) as
\begin{equation} {\rm d}F_\mu \propto 4\mu _d\, \dfrac {\dot {D}}{D}\, \dfrac {1}{r}. \end{equation}
Integrating (4.9) over the entire volume of the pancake, we get an estimate for the total viscous dissipation force working against the internal flow as
\begin{equation} F_\mu \propto \dfrac {4}{3} \pi \mu _d D_0^3 \dfrac {\dot {D}}{D^2}. \end{equation}
A drop under impulsive acceleration starts at zero velocity and then asymptotically accelerates towards a maximum velocity equal to that of the ambient flow (given that the drop is still intact). This acceleration is driven by the pressure and viscous stresses applied by the outside flow on the drop interface, the corresponding magnitude of which is dictated by the instantaneous Reynolds number and relative velocity with the ambient medium
$V_{\mathit{rel}}$
. As the drop continues to accelerate, its relative velocity with respect to the ambient medium continues to decrease, which reduces both the pressure and the viscous stresses. For water--air systems where
$\mathit{Oh}_o$
is low and
$\rho$
is high, the drop does not show significant accelerations and
$V_{\mathit{rel}}$
remains close to the initial relative velocity of
$V_0$
for almost the entire duration of the breakup process. For such cases, setting
$V_{\mathit{rel}}$
equal to
$V_0$
would be a valid assumption. Conversely, if the drop shows significant accelerations and gains velocities that are a substantial fraction of
$V_0$
(which is the case for low
$\rho$
or high
$\mathit{Oh}_o$
systems),
$V_{\mathit{rel}}$
cannot be assumed to be equal to
$V_0$
anymore. It then becomes essential to derive a scaling for
$V_{\mathit{rel}}$
that is valid for the whole parameter space. For this purpose, we utilise the drag equation for a sphere, i.e.
\begin{equation} F_{\mathit{drag}} = \dfrac {1}{2} C_D \rho _o V_{\mathit{rel}}^2 A \end{equation}
to evaluate the drag forces acting on the drop at time
$t=0$
(i.e. when the drop is spherical and at rest):
\begin{equation} \rho _d \frac {\pi }{6} D^3 a_0 = F_{\mathit{drag}}(t=0) = \frac {\pi }{8} C_{D0} \rho _o V_0^2 D^2. \end{equation}
Equation (4.12) gives us a scale for the acceleration experienced by the drop at
$t=0$
:
\begin{equation} a_0 = \frac {3}{4} \frac {1}{\rho D} C_{D0} V_0^2. \end{equation}
A scale describing the velocity of the drop relative to the ambient medium after a time of the order of the deformation time scale
$\tau$
(§ 4.2) can be obtained by evaluating
$V_0 - \tau a_0$
, which when simplified gives
\begin{equation} V_{\mathit{rel}} \propto K_v V_0 \, , \text{ where } K_v = \left ( 1- \frac {3}{4} \frac {C_{D0}}{\sqrt {\rho }} \right )\!. \end{equation}
In figures 5 and 9 we observe that the drop’s velocity decreases with an increase in
$\rho$
and a decrease in
$\mathit{Oh}_o$
. These observations are in line with the relative velocity scale obtained in (4.14), where
$C_{D0}$
is a function of
$\mathit{Re}_0 = \sqrt {\mathit{We}_0}/\mathit{Oh}_o$
. Any increase in the ratio between
$C_{D0}$
and
$\sqrt {\rho }$
leads to a decrease in
$V_{\mathit{rel}}$
. It is important to note that the use of a drag coefficient for a sphere in the derivation of
$V_{\mathit{rel}}$
is a simplification that allows us to obtain an explicit estimate for
$V_{\mathit{rel}}$
. However, this choice results in an underestimation of the average acceleration experienced by the drop, since a pancake shape is less aerodynamic and has a larger frontal area compared with a sphere. This results in an overestimation of
$V_{rel}$
, and hence, an overestimation of pressure forces as well as viscous dissipation. In addition, as the drop accelerates and deforms, its instantaneous Reynolds number also changes, which would affect the drag coefficient. However, the threshold is expected to be overestimated as a result of this simplification, since the aerodynamic forces are proportional to the square of the relative velocity.
Here
$\dot {D}$
is estimated to be equal to
$2V_{\mathit{rel}}/(\sqrt {\rho }+1)$
, described as the Dimotakis velocity (Dimotakis Reference Dimotakis1986) and also used by Marcotte & Zaleski (Reference Marcotte and Zaleski2019) in their work. Finally, in all simulations conducted in the current work across the parameter space, it was observed that the aspect ratio
$A_{xr}$
at
$t^*\approx 1$
was almost always
$0.175$
, and we use this value and the corresponding
$D$
for the calculation of the forcing scales.
At this stage, we have derived all the relevant forcing scales and related parameters for deriving a new parameter describing the fragmentation threshold. However, this derivation hinges upon an important assumption that the pancake is a flat disk for all the cases. This is however not true for all the high
$\mathit{Oh}_o$
cases that form a forward-facing pancake (concave shape pointing downstream). This difference in geometry affects the potential flow and the corresponding pressure field around the drop, the pressure jumps due to surface curvature and the corresponding internal flow vectors.
It is also important to note that the use of the deformation time scale
$\tau$
to estimate the velocity scale implicitly inherits the assumption that the drop flattens into a pancake in time
$\tau$
, and the balance of forces acting on the drop control whether the pancake invariably develops a bag. Hence, we expect
$C_{\mathit{breakup}}$
to correctly capture the threshold only for drops that go through a critical pancake shape during their deformation process.
Equation (4.3) can now be rewritten incorporating the pressure deficit and stagnation pressure forces (4.7), surface tension force (4.4) and viscous resistance in a drop fluid (4.10), i.e.
\begin{equation} C_{\mathit{breakup}} = \frac { F_{p,d} }{ F_{p,s} + F_\sigma + F_\mu }, \end{equation}
which finally results in
\begin{equation} C_{\mathit{breakup}} = \dfrac { (a^2/48)\,\mathit{We}_0 K_v^2 D^2 }{ 2D^2(1+A_{xr}) + (4/3)\mathit{Oh}_d\, (\rho \mathit{We}_0)^{0.5}\dot {D} + (1/3)\mathit{We}_0 K_v^2 }. \end{equation}
For a water--air analogous systems,
$\mathit{Oh}_d\lt \lt 1$
and
$\rho \gt 500$
, which implies that
$K_v \approx 1$
(
$C_{D0}\approx 1$
for large
$\mathit{Re}_0$
) and
$\dot {D} \approx V_{\mathit{rel}}/\sqrt {\rho }$
. Then
$C_{\mathit{breakup}}$
can be simplified to
\begin{equation} C_{\mathit{breakup}} = \dfrac { (a^2/48)\,\mathit{We}_0 D^2 }{ 2D^2(1+A_{xr}) + (1/3)\mathit{We}_0 }. \end{equation}
For obtaining
$C_{\mathit{breakup}}$
using (4.16), an explicit equation for the drag coefficient is required. We use the following relationship provided by Turton & Levenspiel (Reference Turton and Levenspiel1986):
\begin{equation} C_{D0} = \frac {24}{\mathit{Re}_0} \big( 1 + 0.173 \mathit{Re}_0^{0.657} \big ) + \frac {0.413}{1+16300 \mathit{Re}_0^{-1.09}}. \end{equation}
Note that
$C_{D0}$
is a function of
$\mathit{Re}_0$
but can also be expressed in terms of
$\mathit{Oh}_o$
and
$\mathit{We}_0$
since
$\mathit{Re}_0 = \sqrt {\mathit{We}_0}/\mathit{Oh}_o$
.
The
$\mathit{We}_{cr}$
vs
$\mathit{Oh}_d$
plot in figure 15 is recreated in figure 18(c) using
$C_{\mathit{breakup}}$
(4.16) on the
$y$
axis. When scaled according to this new non-dimensional parameter, almost all simulation points move to a narrow range of
$0.083\lt C_{\mathit{breakup}}\lt 0.137$
with a standard deviation of
$0.02$
about a mean of
$0.11$
, when compared with
$10\lt \mathit{We}_{cr}\lt 74$
in the former plot. The highest
$\mathit{We}_{cr}$
cases corresponding to
$\mathit{Oh}_d$
,
$\mathit{Oh}_o=0.1$
and
$\rho =10,50$
in figure 15 when plotted according to (4.16) achieve substantially lower threshold values, much more in line with the
$C_{\mathit{breakup}}$
values obtained for other simulation cases. It is also observed that the non-trivial fragmentation morphologies such as backward bag-plume and backward-transition breakup cases show relatively the largest
$C_{\mathit{breakup}}$
values among all other cases with the same
$\mathit{Oh}_d$
.
We also plot
$C_{\mathit{breakup}}$
values of all the available experimental data for critical backward bag breakup as a set of reference data for the plot. As expected, the experimental data remains bound within a narrow extent of
$C_{\mathit{breakup}}$
values, similar to its
$\mathit{We}_{cr}$
analogue. The
$C_{\mathit{breakup}}$
values do start to drop off as we approach larger and larger
$\mathit{Oh}_d$
values, i.e.
$C_{\mathit{breakup}}$
appears to overestimate the effort required to achieve critical breakup in very high
$\mathit{Oh}_d$
drops compared with experiments. This is in line with the observations in previous experimental works (Hinze Reference Hinze1955; Hsiang & Faeth Reference Hsiang and Faeth1995) where it was observed that bag breakup becomes progressively difficult for higher drop Ohnesorge numbers and ultimately stops manifesting for
$\mathit{Oh}_d\gt 2$
. Other breakup modes such as multimode and sheet thinning become the critical breakup modes for very high
$\mathit{Oh}_d$
systems. The inherent assumption in the derivation of
$C_{\mathit{breakup}}$
was to assume that critical breakup morphology is a bag breakup, i.e. a drop first flattens into a pancake that then blows in to a bag. Any major deviation from this general breakup mechanism is expected to drastically change the deformation and breakup physics of such drops. Hence, as
$\mathit{Oh}_d$
increases past
$1$
, physics related to multimode and sheet thinning breaks starts to become significant and, hence, must be considered in the derivation of any parameter attempting to define the critical breakup criteria.
Another factor worth considering is the initial Reynolds number required (i.e.
$\mathit{Re}_0 = \sqrt {\mathit{We}_0}/\mathit{Oh}_o$
) for very large
$\mathit{We}_0$
values to show a breakup for high
$\mathit{Oh}_d$
drops. For very large
$\mathit{We}_0$
values, very large
$\mathit{Re}_0$
values are expected, which would make the external flow more chaotic, which can have a major impact on the pressure forces experienced by the drop, its interaction with downstream vortices and the intensity of surface waves (§ 3.2). These effects have not been taken into account in our derivation of
$C_{\mathit{breakup}}$
.
For the non-dimensional parameter space considered in the current work,
$C_{\mathit{breakup}}$
captures the dominant physics very well and succeeds in compressing the rather large variance in
$\mathit{We}_{cr}$
values observed due to very low
$\rho$
and high
$\mathit{Oh}_d$
and
$\mathit{Oh}_o$
values. Therefore, for the purposes of the work,
$C_{\mathit{breakup}}$
fulfils our requirements.
5. Conclusions
The current work aimed to clarify two major questions: (i) the effect of
$\mathit{We}_0$
,
$\rho$
,
$\mathit{Oh}_o$
and
$\mathit{Oh}_d$
on drop deformation and breakup characteristics, at or near the
$\mathit{We}_{cr}$
and for physical properties that are vastly different from commonly available water--air systems; and (ii) the effect of these non-dimensional parameters on drop critical breakup morphology.
As explained extensively in § 1, most of the currently accepted ideas on threshold secondary atomisation, such as the independence of
$\mathit{We}_{cr}$
with respect to
$\mathit{Oh}_d$
given it is below
$0.1$
, or the critical breakup morphology being solely a backward bag, etc. originate from experimental works done for a small parametric space occupied by water--air analogous systems. The current work aimed to shed some light on these accepted ideas and provide a more complete picture of the process of secondary atomisation of Newtonian drops. For this purpose, a parametric sweep across all the involved non-dimensional parameters, i.e.
$\{\rho , \mathit{Oh}_o, \mathit{Oh}_d\}$
(§ 2.4) was performed using axisymmetric simulations on Basilisk, with
$\mathit{We}_0$
values increased until a non-vibrational (critical) breakup was achieved for a given
$\{\rho , \mathit{Oh}_o, \mathit{Oh}_d\}$
set. The large amount of simulation data obtained from the comprehensive parameter exploration resulted in the following key findings.
-
(i) The internal flow inside a drop away from its core can be motivated by two forces: the pressure difference between its poles (
${\Delta P}_{\textit{dri}v\textit{e}}$
) and its periphery, and the shear stresses acting on its upstream surface. The competition between the two and the relative strengths of the two dictates whether the the drop pancake is flat or forward-facing. An internal flow predominantly driven by a pressure difference results in a flat pancake, whereas a shear-stress driven internal flow results in a forward-facing pancake (§ 3). -
(ii) Density ratio controls
${\Delta P}_{\textit{dri}v\textit{e}}$
, whereas
$\mathit{Oh}_o$
controls the shear stresses on the upstream surface. Thus, both these non-dimensional parameters are important when predicting the orientation of the pancake. -
(iii) Given an external forcing, the sensitivity (i.e. local accelerations) of the different parts of the drop when experiencing external forces is directly controlled by the local inertia of that region. The morphology of the drop (forward or backward bag) as it continues to deform past the pancake stage is then controlled by the relative accelerations of different parts, which depends on the dominant mechanism driving the internal flow (pressure versus shear stresses), the rate of evacuation of the core, the strength of downstream vortices and the density ratio of the drop fluid.
-
(iv)
$\mathit{Oh}_o$
controls the length scales and time scales of the vortices produced in the wake of the drop. If the vortices do not detach from the periphery of the drop, higher local accelerations at the periphery are observed, which results in the formation of a forward-facing bag. This is generally only the case for low
$\rho$
drops. -
(v) Given
$\rho$
and
$\mathit{Oh}_o$
are the same, the drop Ohnesorge number
$\mathit{Oh}_d$
controls the response of the drop fluid to the external forcing. A low
$\mathit{Oh}_d$
may allow the development of a surface instability at the upstream pole, resulting in the formation of a plume, and thus making a backward bag-plume morphology as the threshold morphology. The drop Ohnesorge number
$\mathit{Oh}_d$
also controls the rate of evacuation of the core, with lower
$\mathit{Oh}_d$
cases showing slower evacuation and, thus, motivating the formation of a plume. -
(vi) The critical Weber number
$\mathit{We}_{cr}$
over the explored parameter space is obtained and plotted against
$\mathit{Oh}_d$
(figure 15) to recreate the plot presented in Hsiang & Faeth (Reference Hsiang and Faeth1995). It is found that
$\mathit{We}_{cr}$
varies significantly in the space
$\mathit{Oh}_d\lt 0.1$
, achieving values as high as
$70$
for the lowest density ratio and ambient Ohnesorge number systems. Furthermore, backward bag-plume and forward bag morphologies are observed for threshold fragmentation in addition to the trivial backward bag morphology. -
(vii) Based on all the insights gained from the simulation results, a phase diagram (figure 16) is constructed describing the various deformation pathways a spherical drop may undertake under impulsive acceleration depending on the properties. The drop deformation path shows its first split at the shape of the pancake, which then leads to the possible threshold fragmentation morphologies.
-
(viii) The simulations also allow us to explore bag inflation characteristics for backward bag breakups for a greater parameter space, and extract some general conclusions on the influence of each parameter in modulating the growth rate and the associated time scales. Additionally,
$t^* \approx 1$
is found to be a good measure for the time scale required for the initiation of bag inflation. -
(ix) Finally, a non-dimensional parameter (
$C_{\mathit{breakup}}$
) is derived by finding the ratio of forces supporting the fragmentation of the drop to the forces opposing it. Surface tension forces, pressure forces and viscous dissipation are considered. The change in velocity of the drop relative to the ambient medium (which can be significant for low
$\rho$
drops) is accounted for when calculating the forces. The term
$C_{\mathit{breakup}}$
thus obtained provides a more complete alternative to the Weber number as an indicator of drop threshold criteria in the dimensional space of the current study, by capturing the effects of all three studied non-dimensional parameters, i.e.
$\rho$
,
$\mathit{Oh}_o$
and
$\mathit{Oh}_d$
, on drop deformation and breakup.
Acknowledgement
Computational resources provided by the COVID-19 HPC Consortium through time on Blue Waters (NCSA) were used for the simulations.
Funding
S.D. and A.P’.s participation was supported by the DOE Office of Science through the National Virtual Biotechnology Laboratory (NVBL), a consortium of DOE national laboratories focused on response to COVID-19, with funding provided by the Coronavirus CARES Act. We would also like to thank Texas Advanced Computing Center for providing computing time on the NSF funded Frontera supercomputer, through the Pathways grant program. The computing time allowed us to simulate the 3-D cases that were used to validate the accuracy of the axisymmetric simulations.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical scheme
Basilisk solves incompressible Navier--Stokes multiphase flow equations ((A1), (A2) and (A3)) on a quad/octree discretised grid, which allows variable mesh densities at the interface (Popinet Reference Popinet2003) and, therefore, accurately captures capillary scale phenomena. Conservation of momentum and mass results in the equations
\begin{align} \rho (\partial _t \boldsymbol{u} + \boldsymbol{u}.\boldsymbol{\nabla u}) &= -\boldsymbol{\nabla } p + \boldsymbol{\nabla }.(2 \mu \boldsymbol{D}) + \sigma \kappa \delta _s \boldsymbol{n}, \end{align}
\begin{align} \partial _t \rho + \boldsymbol{\nabla }.(\rho \,\boldsymbol{u}) &= 0, \end{align}
\begin{align} \boldsymbol{\nabla .u} &= 0, \end{align}
where
$\boldsymbol{u} = (u, v, w)$
is the fluid velocity,
$\rho \equiv \rho (\boldsymbol{x}, t)$
is the fluid density,
$\mu \equiv \mu (\boldsymbol{x}, t)$
is the dynamic viscosity and
$\boldsymbol{D}$
is the deformation tensor defined as
$D_{ij} = (\partial _i u_j + \partial _j u_i)/2$
. The Dirac distribution function
$\delta _s$
allows inclusion of surface tension forces in the momentum equation by switching on the surface tension term only at the interface between the fluids;
$\sigma$
is the surface tension coefficient,
$\kappa$
and
$\boldsymbol{n}$
the curvature and the normal to the interface, respectively. We compute
$\kappa$
using the height function formulation as described by Torrey et al. (Reference Torrey, Cloutman, Mjolsness and Hirt1985), with attention given to address under-resolved interfaces. The surface tension term is calculated using the continuum surface force approach first described in Brackbill, Kothe & Zemach (Reference Brackbill, Kothe and Zemach1992), with special care taken to ensure that the conditions described in Francois et al. (Reference Francois, Cummins, Dendy, Kothe, Sicilian and Williams2006) are satisfied to prevent parasitic currents.
To maintain the single equation formulation of the momentum equation, the two fluids are represented using a volume fraction
$c(\boldsymbol{x}, t)$
according to which
$\rho$
and
$\mu$
are defined as
\begin{align} \rho &= c \rho _1 + (1-c) \rho _2, \end{align}
\begin{align} \mu &= c \mu _1 + (1-c)\mu _2, \end{align}
where
$\rho _1, \rho _2$
and
$\mu _1, \mu _2$
are the densities of the first and second fluid in the domain, respectively. In this formulation, the density advection equation is replaced with a volume fraction advection equation:
\begin{align} \partial _t c + \boldsymbol{\nabla } \boldsymbol{\cdot}(c\,\boldsymbol{u}) = 0. \end{align}
The entire computational domain is discretised using squares for two dimensions (quadtree) and cubes for three dimensions (octree) and then organised in a hierarchy of cells. The mesh resolution is adaptive in nature, and hence, the two-fluid interface is resolved at a much higher resolution than other computationally less interesting regions of the domain. This allows for large savings in the computational costs for two-phase simulations. Any cell serving as a parent computational element can undergo further refinement into four or eight equal children cells for two- and three-dimensional computations, respectively. Each of the children cells, in turn, can act as a parent cell if further refinement is warranted. This successive refinement continues until a (user-defined) threshold criterion for error is satisfied or a maximum refinement level is reached. A wavelet-based error estimation is used to estimate errors associated with the specified fields (Popinet Reference Popinet2015; van Hooft et al. Reference van Hooft, Popinet, van Heerwaarden, van der Linden, de Roode and van de Wiel2018). The maximum allowed refinement, corresponding to the smallest allowed cell size, is constrained by a user-specified minimum allowed cell dimension, which is defined by a parameter called ‘maximum level’, a maximum level of
$N$
corresponding to a minimum cell size of
$L/2^N$
.
Appendix B. Choice of numerical parameters
Before using Basilisk for the production runs, it is essential to test for independence of the solutions with respect to both the maximum mesh resolution (normally achieved at the interface), wavelet-error thresholds for relevant field variables and the tolerance of the Poisson solver. For drop simulations, the accuracy of the calculated interface and the velocity fields must be ensured for correct retrieval of surface stresses, and correspondingly the temporal development of the drop deformation process. Hence, in the current simulations the extent of adaptive mesh refinement performed at a computational cell is restricted by the allowable maximum wavelet errors for velocity (
$\chi _u$
) and volume fraction (
$\chi _c$
) fields. Additionally, a maximum allowed refinement level (
$N$
) is specified that enforces a strict allowed minimum cell size of
$L/2^N$
across all the computational cells. This helps prevent unbounded mesh refinement if the error threshold criteria converges slowly with a decrease in cell size in a computational cell, which in turn eases the restriction on the simulation time step (dependent on the cell with the highest Courant--Friedrichs--Lewy number in the domain). Finally, the tolerance of the Poisson solver (
$\epsilon$
) is also a critical parameter that can affect the accuracy of the velocity field and, hence, the interface.
We test the independence of the axisymmetric simulations with respect to
$\chi _u$
,
$N$
and
$\epsilon$
using a test case with
$\{\rho , \mathit{Oh}_o, \mathit{Oh}_d, \mathit{We}_0\} = \{500, 0.001, 0.01, 10\}$
on an axisymmetric domain (figure 2) with
$L=16$
,
$D=1$
and
$V_0=1$
. A drop-ambient system with these physical properties is expected to deform close to fragmentation, but instead of fragmenting it retracts towards its neutral spherical shape. Thus, it provides a good representative case for identifying any spurious diffusion accumulated in the system. For testing the influence of
$\chi _c$
, we simulate a case with
$\rho , \mathit{Oh}_o, \mathit{Oh}_d, \mathit{We}_0 = 500, 0.01, 0.1, 16$
, which is expected to inflate into a bag morphology and fragment, and hence, is a good representative case for the accuracy of the interface calculation. The default values of
$\chi _c$
,
$\chi _u$
,
$\epsilon$
and
$N$
for all cases is set to
$10^{-6}$
,
$10^{-3}$
,
$10^{-4}$
and
$13$
, respectively, and the parameter of concern is varied to test the independence of the solution with respect to it.
Volume-of-fluid plots at specific times for different error thresholds for the volume fraction field
$c$
(
$\chi _c$
) and velocity fields (
$\chi _u$
). The role of the tolerance of the Poisson solver (
$\epsilon$
) is also shown through its effect on the interface.
The
$x$
coordinates of the centre of mass, the
$x$
velocity of the centre-of-mass and the axis ratio are shown for different thresholds for wavelet errors of the volume fraction fields (
$\chi _c$
) (a--c), different threshold for wavelet error of the velocity field (
$\chi _u$
) (d--f), and different tolerances of the Poisson solver (g--i), respectively. In (a--c) the maximum allowed refinement level
$N$
refers to a minimum cell size of
$L/2^N$
. Thus, 256, 512 and 1024 cells per diameter correspond to
$N=12$
,
$N=13$
and
$N=14$
respectively, given
$L=16$
and
$D=1$
.
In figure 19 the drop interface at specific times is shown for different values of
$\chi _c$
,
$N$
,
$\chi _u$
and
$\epsilon$
. Here
$\chi _c$
is varied across
$10^{-3}$
,
$10^{-6}$
and
$10^{-9}$
, and we observe that the interface is identical for all three values. This indicates that the interface is converged with respect to
$\chi _c$
. The effect of
$N$
however is significant, with the drop showing far more numerical diffusion for
$N=12$
compared with
$N=13$
and
$N=14$
. This results in the drop achieving less deformation and, thus, starts retracting earlier for
$N=12$
. The effect of
$\chi _u$
and
$\epsilon$
is negligible for the two time instances shown, with the interface being identical for all the explored values of
$\chi _u$
and
$\epsilon$
.
Figure 20 shows the effect of
$N$
,
$\chi _u$
and
$\epsilon$
on the
$x$
coordinate of the centre of mass,
$x$
velocity of the centre of mass and the aspect ratio of the drop, resulting in three sets of plots for each parameter. Plots for
$\chi _c$
are not shown since we have already established that the interface is converged with respect to
$\chi _c$
. It is also found that the computational costs associated with the three
$\chi _c$
values are within
$1\,\%$
of each other. This makes the choice of
$\chi _c$
less critical, and hence, the default value of
$\chi _c=10^{-6}$
is chosen for the production runs.
Figure 20(a,b,c) shows the effect of
$N$
on the
$x$
coordinates of the centre of mass,
$x$
velocities of the centre of mass and aspect ratios of the drop, respectively. It is observed that the three
$N$
values start to diverge as the drop reaches its minimum aspect ratio, with
$N=12$
showing the highest values corresponding to the least deformation. As the maximum allowed refinement increases, the drop starts to deform more, with
$N=14$
showing the highest deformation. This points to the larger numerical diffusion in the
$N=12$
case. The
$N=13$
and
$N=14$
cases show much closer results for all three plots. We find that
$N$
also has a dramatic effect on computational costs, with
$N=14$
requiring approximately three times the computational time as required for
$N=13$
. The
$N=13$
case, which is equivalent to
$512$
cells per diameter for
$L=16$
, provides a good balance between accuracy and computational costs. For all simulations in this study,
$N=13$
is chosen as the default value for
$L=16$
. For the highest
${Re}_0$
cases (
$\mathit{Oh}_o=0.0001$
),
$N=14$
is chosen to ensure that the low viscosity ambient flow is resolved accurately.
Figure 20(d,e,f) plots
$x_{cm}$
,
$V_{cm,x}$
and
$A_{xr}$
for four different values of
$\chi _u$
. We observe that apart from the lowest value of
$\chi _u=10^{-2}$
, the three properties are almost identical for all values of
$\chi _u$
. The
$\chi _u=10^{-2}$
case shows a lower
$A_{xr}$
, resulting from higher dissipation of the energy contained in the drop, compared with the other three cases, resulting in a slightly higher
$x_{cm}$
and
$V_{cm,x}$
. The results for
$\chi _u=10^{-4}$
and
$\chi _u=10^{-5}$
are almost identical, with the computational cost for
$\chi _u=10^{-4}$
being approximately 2.5 times less than that for
$\chi _u=10^{-5}$
. Hence,
$\chi _u = 10^{-4}$
is chosen for the production runs.
Finally, from the convergence plot for the tolerance of the Poisson solver in figure 20(g,h,i), we observe that the
$x$
coordinates of the centre of mass,
$x$
velocities of the centre of mass and axis ratios are almost identical for all three values of
$\epsilon$
. We hence use the lowest value of
$\epsilon =10^{-3}$
for the production runs, as it provides a converged result with the least computational cost.
Appendix C. Comparison to 3-D simulations
The fluid interface for different time instances for axisymmetric (shown in black) and 3-D (shown in red) simulations for five different cases, with the drop-ambient system properties mentioned in the figure as
$\{\rho , \mathit{Oh}_o, \mathit{Oh}_d, \mathit{We}_0\}$
. Two of these cases, shown in (a) and (b), are chosen to verify if the axisymmetric simulations can capture the formation of a forward pancake in two different contexts: a large density-ratio drop (
$\rho =500$
) and a low density-ratio drop (
$\rho =10$
). The cases shown in (c) and (d) help us verify the physical validity of forward bags observed during threshold fragmentation of low
$\rho$
drops and small
$\mathit{Oh}_o$
values. The final case (e) serves to highlight the differences between axisymmetric and 3-D simulations for cases where both the drop and ambient Ohnesorge numbers are the smallest. The pressure fields at
$t^*\approx 1$
are also shown.
To justify the use of axisymmetric simulations for the wide parameter space of density ratios and viscosities under consideration, we conduct 3-D simulations for some benchmark cases in order to verify if 3-D simulations produce similar non-trivial morphologies as the axisymmetric simulations. A qualitative match with axisymmetric simulations would justify their use for prediction in the respective pancake and bag orientations for the parameter space. Jain et al. (Reference Jain, Tyagi, Prakash, Ravikrishna and Tomar2019) states that, ‘For the drops with high
$\rho$
, flow around the drop has relatively low effect on the drop deformation, morphology and the breakup’. If this is the case, the ideal benchmark drop-ambient systems should have either low
$\rho$
drops in low
$\mathit{Oh}_o$
flows, or high
$\rho$
drops in high
$\mathit{Oh}_o$
flows in order to magnify the role of ambient flow on the drop. This is a computationally favourable scenario for us, since both low
$\rho$
and low
$\mathit{Re}_0$
systems have significantly reduced computational requirements. Keeping this in mind, we specify the following six benchmark cases in three dimensions.
-
(i) Two of the cases are designed to verify the formation of forward pancakes in axisymmetric simulations, one is a large density-ratio drop (
$\rho =1000$
) in a shear-stress dominated (large
$\mathit{Oh}_o$
) system, another is a low density-ratio drop (
$\rho =10$
) in a pressure dominated system. These simulations are run at the same refinement level as the axisymmetric simulations. -
(ii) Two of the cases test the ability of axisymmetric simulations to correctly predict the transition from a forward pancake to forward bag for low density-ratio drops in high
$\mathit{Re}_0$
systems. Both the simulations are run at the same refinement level as the axisymmetric simulations. -
(iii) An additional case at an intermediate density ratio of
$\rho =100$
is used to highlight the differences between the two simulation types for large ambient Reynolds number cases (
$\mathit{Oh}_o=0.0001$
) and a low viscosity drop (
$\mathit{Oh}_d=0.001$
). This choice renders the ambient flow non-axisymmetric since the formation of turbulent vortices is a purely 3-D phenomenon, making such a flow system impossible to perfectly reproduce with only axisymmetric simulations. The high
$\mathit{Re}_0$
of the flow coupled with an intermediate density ratio makes the simulations computationally extremely expensive making a 3-D simulation until fragmentation unfeasible. These simulations have been run only until
$t^* \approx 1$
, i.e. until the formation of the pancake.
Backward bag formation has not been tested here, since it has already been considered in § 2.3. The fluid interface for axisymmetric and 3-D simulations for the cases described in (i) and (ii) are shown in figure 21(a–d), with a red colour representing the 3-D simulations. All four cases show a good qualitative match between the two simulation types, with the formation of a forward pancake observed in all cases. The cases with
$\mathit{Oh}_o \leqslant 0.001$
(shown in panels a and b) show an excellent match in the interface shape all throughout the deformation process. For the
$\mathit{Oh}_o=0.0001$
cases (shown in panels cand d), the interface shape differs more, with the 3-D simulations showing a more pronounced forward pancake and a larger forward bag compared with the axisymmetric simulations. This can be attributed to the differences in viscous dissipation between axisymmetric and 3-D simulations. However, most importantly, the 3-D simulations transition from a forward pancake to a forward bag in a manner similar to axisymmetric simulations, thus supporting the physical validity of forward bag morphology observed in the axisymmetric simulations.
Figure 21(e) shows the pressure field at
$t^* \approx 1$
for the
$\mathit{Oh}_o=0.0001$
case for
$\rho =100$
, as described in (iii). It is observed that the pressure fields even for such high
$\mathit{Re}_0$
cases are close to toroidal and interact similarly with the periphery of the drop. However, the downstream stagnation pressures are much higher for the axisymmetric simulations. This is consistent with the observations of Ling & Mahmood (Reference Ling and Mahmood2023). Thus, for all very large
$\mathit{Re}_0$
simulations with low
$\mathit{Oh}_d$
performed in this work, both the pancake shape and subsequent fragmentation morphology are not well justified. All simulation results for such cases are to be considered with this caveat in mind. It is hypothesised that highly turbulent ambient flow in low
$\mathit{Oh}_o$
systems produces non-axisymmetric interfacial perturbations that are not reproduced in the axisymmetric simulations. These perturbations however have little significance for the high
$\mathit{Oh}_d$
cases, since they are rapidly dissipated at time scales smaller than the deformation time scale. However, for low
$\mathit{Oh}_d$
drops, these surface perturbations persist for longer time scales and, hence, can significantly influence the subsequent deformation. Thus, the axisymmetric simulations are not expected to accurately predict the fragmentation morphology for cases with high ambient Reynolds numbers coupled with very low drop viscosities.
Open access
