Density contrast matters for drop fragmentation thresholds at low Ohnesorge number

TL;DR

This study numerically investigates how the density contrast between liquid and gas influences drop fragmentation thresholds at low Ohnesorge numbers, revealing that density contrast significantly affects the critical Weber number for breakup regimes.

Contribution

It introduces a theoretical prediction for the transitional Weber number based on density contrast, highlighting the stabilizing effect of small density contrasts in drop fragmentation.

Findings

Density contrast modifies the critical Weber number for breakup.

Small density contrasts have a stabilizing effect explained by inertial effects.

Theoretical model predicts transition Weber number as a function of density contrast.

Abstract

We address numerically the deformation and fragmentation dynamics of a single liquid drop subject to impulsive acceleration by a unidirectional gas stream. The density ratio between the liquid and gaseous phases is varied in the $10-2000$ range and comparisons are made with recent drop breakup experiments. We show for low Ohnesorge numbers that the liquid-gas density contrast significantly modifies the critical Weber number for the transition between bursting and stripping fragmentation regimes, on the one hand, and for drop fragmentation, on the other hand. We suggest a simple theoretical argument to predict the transitional Weber number as a function of the density contrast and show that the stabilising influence of small contrasts can be explained by inertial effects in the nonlinear coupling between drop stretching and centroid acceleration.

Equations

$${\rho }_L\mathcal{V}{U}_{\ast }={\rho }_L\mathcal{V}{U}_1+{\rho }_G{C}_M\mathcal{V}{U}_1,$$

$$U_1=U_{\ast }{\left(1+\frac{1}{2}{\mathtt{r}}_{\mathtt{d}}^{-1}\right)}^{-1},$$

$${\rho }_L{\dot{R}}^2\sim \sigma \kappa ,$$

$$\frac{{\rho }_{L}e{\dot{R}}^2}{2\sigma }\sim 1,$$

$$1\sim \text{We}\frac{e}{2R_0}\sim \text{We}\frac{2}{3}\frac{R_0^2}{R^2},$$

$$1\sim \text{We}\frac{2}{3}(1+t_{\times }^{\ast }{)}^{-2},$$

$${\mu }_G\frac{U_0}{\delta }\sim {\mu }_L\frac{U_L}{e},$$

$$t_S^{\ast }={{\mathtt{r}}_{\mathtt{d}}}^{-\frac{1}{2}}\frac{U_0}{R_0}{t}_S\sim {{\mathtt{r}}_{\mathtt{d}}}^{-\frac{1}{2}}{\mathtt{r}}_{\mathtt{v}}\sqrt{\frac{1}{\text{Re}}}.$$

$$\text{We}\sim \frac{3}{2}{\left(1+\lambda {{\mathtt{r}}_{\mathtt{d}}}^{-\frac{1}{2}}{\mathtt{r}}_{\mathtt{v}}\sqrt{\frac{1}{\text{Re}}}\right)}^2,$$

$${\rho }_L\frac{d}{dt}\left(U_C\mathcal{V}\right)=C_d{\rho }_G(U_0-U_C{)}^2\pi R^2,$$

$$r{\partial }_{t}e+{\partial }_r(r{\bar{u}}_r)=0,$$

$$p_L(R)-p_L(0)\sim -\frac{{\alpha }^2}{8}{\rho }_G(U_0-U_C{)}^2{R}^2/R_0^2+\frac{3}{2}\sigma R^2/R_0^3.$$

$${\rho }_L({\partial }_t{u}_r+u\cdot \nabla u_r)=-{\partial }_{r}p.$$

$$\frac{\ddot{R}}{R}\ \sim\,{\tt r_{d}}^{-1}\left(\frac{\alpha^{2}}{4}\,\Big{(}1-\frac{U_{C}}{U_{0}}\Big{)}^{2}-\frac{3}{\text{We}}\right)\,\frac{U_{0}^{2}}{R_{0}^{2}},$$

$$\dot{U}_{C}\ \sim\ \frac{3}{4}C_{d}\,{\tt r_{d}}^{-1}\Big{(}1-\frac{U_{C}}{U_{0}}\Big{)}^{2}\,\frac{U_{0}^{2}}{R_{0}^{3}}\,R^{2},$$

Full text

Drop fragmentation is the complex phenomenon by which single liquid drops break into child droplets of varying sizes, occurring in a large diversity of natural events (e.g. rain) or industrial processes (e.g. fuel injection in propulsion systems). It is also referred to as secondary atomisation in the context of the classical jet atomisation problem where, through the development of successive hydrodynamic instabilities, a high-speed liquid jet destabilises into thin liquid sheets prone to fingering and (primary) drop ejection. The further fragmentation of the resulting drops, if any, depends on the flow conditions and eventually determines the droplets size distribution, a subject of many active debates and a key element toward control of a spray mixing rate Eggers08 ; VB09 ; ling16 .

The fundamental mechanisms underlying drop fragmentation have been addressed by a vast amount of experimental and numerical studies, which considered the breakup of a single, initially spherical drop immersed in a unidirectional high-speed gas stream. Drop breakup experiments carried out in the last few decades are essentially twofold: a liquid drop is either impulsively accelerated by a shock-wave in a dedicated tube Hsiang95 ; Faeth95 ; chou98 ; dai01 or released sideways in a high-speed cross-flow liureitz97 ; cao07 ; zhao2010 ; opfer14 . These experiments have established the existence of a variety of deformation and breakup regimes depending on the flow conditions Hinze55 ; ranger69 ; Gelfand74 ; krzeczkowski80 ; pilch87 ; wierzba90 . Due to practical constraints however, drop breakup experiments have mostly focused on liquid drops whose density $\rho_{L}$ was large compared to that of the ambient gas $\rho_{G}$, with density contrasts ${\tt r_{d}}=\rho_{L}/\rho_{G}$ typically larger than 500 (as in the case of a water-air or ethanol-air system at ambient temperature and atmospheric pressure opfer14 ). Nevertheless, small density contrasts (${\tt r_{d}}\leq 10$) are also relevant for combustion sprays due to the high fuel injection pressures involved e.g. in rocket propulsion systems Heywood88 . Moreover, information such as the thickness of the sheet of inflated drops remains elusive as yet, due to visualisation difficulties.

Such information, on the other hand, could be easily accessed by means of Direct Numerical Simulations (DNS). Yet numerical studies of drop breakup are limited by the considerable diversity of spatial scales involved as the drop undergoes dramatic stretching and bursting, resulting in huge computational costs, and have long been restricted to small-to-moderate density contrasts with ${\tt r_{d}}$ typically in the $1-80$ range zaleski95 ; han01 ; kekesi14 ; yang16 ; jalaal14 . As a result there is so far little overlap between the flow regimes explored in experiments on the one hand, in simulations on the other hand; and whereas the influence of flow conditions on secondary atomisation has received most attention (especially in terms of the Weber number defined below), this discrepancy arguably contributes to obscure that of the density contrast aalburg ; yang16 . Early simulations zaleski95 ; han01 already suggested that the deformation modes exhibited by the liquid drop at low $\tt r_{d}$ strongly differ from those observed for similar flow conditions at large ${\tt r_{d}}\sim 1000$: therefore it appears desirable to bridge the gap between these two parameter regimes and address the sensitivity of the drop breakup dynamics to the density contrast as the latter is varied over a broad range.

Also of particular interest for industrial applications is to determine the threshold (in terms of flow conditions) beyond which liquid drops become unstable to fragmentation. In the literature this threshold is characterised by means of two dimensionless parameters, namely the Weber number $\text{We}=\rho_{G}R_{0}U_{0}^{2}/\sigma$ quantifying the ratio of inertial to surface tension forces, and the Ohnesorge number quantifying the effet of liquid viscosity ${\text{Oh}}=\mu_{L}/\sqrt{\sigma\rho_{L}R_{0}}$. (Importantly, the latter has been however reported to matter only if exceeding a typical value $\text{Oh}>0.1$.) Here $U_{0}$ is the velocity difference between the gas stream and the drop, $R_{0}$ the (undeformed) drop radius, $\mu_{G}$ the dynamic viscosity of the gas and $\sigma$ the surface tension associated with the gas-liquid interface. Various breakup criteria involving these two parameters are commonly used in non-DNS (Eulerian or Lagrangian) CFD models, where drop breakup is handled with a statistical approach. These numerical criteria either rely on experimental observations (such as the Schmehl model schmehl inferred from pilch87 ; Hsiang95 ) or simple theoretical arguments (such as the Taylor-Analogy Breakup criterion orourke87 , which uses an analogy between the distorted drop and a spring-mass system Taylor63 ). In a fundamental context, VB09 addressed the fragmentation of a falling raindrop and suggested a simple inviscid model to account for the critical Weber number characterising the transition from oscillatory, reversible deformation to exponential stretching and breakup, which they found to be ${\text{We}_{c}}\sim 3$, in excellent agreement with available experiments at large density contrasts.

None of these stability criteria, whether empirical (and based on experiments, hence typically implying ${\tt r_{d}}>500$) or theoretical (and based on simplified models), involves the density contrast between the gaseous and liquid phases. Yet, the important morphological discrepancy between deformation regimes observed at low $\tt r_{d}$ (numerically) and those observed at large $\tt r_{d}$ (mostly experimentally, although simulations at ${\tt r_{d}}=1000$ have been very recently reported for Weber numbers well above the critical value for breakup jain15 ; gaurav2018 ) suggest that the effect of density contrast on the fragmentation threshold remains an open question, which it is the aim of the present paper to investigate.

I Numerical model

I.1 Flow configuration

Simulations were performed using the flow solver Gerris popinet03 ; popinet09 in axisymmetric geometry to model the impulsive acceleration of a spherical liquid drop behind a planar shock-wave. Gerris solves for the two-phase, incompressible Navier-Stokes equations using a finite-volume method with dynamical adaptive mesh refinement, which allows for capturing thin flow structures while mitigating computational costs. The interface is advected using a Volume of Fluid method, and the incompressibility constraint is enforced at each timestep by projecting the velocity field onto the divergence-free manifold using a classical Chorin’s algorithm Chorin67 ; Chorin68 coupled with a multigrid Poisson solver.

Our computational domain is a cylindrical channel with dimensionless radius $R_{cyl}=5$ and length $L_{cyl}=15$ in the streamwise direction, where the unit length is the drop initial radius $R_{0}$ (the channel length was increased whenever necessary up to $L_{cyl}=25$ to observe the successive deformation stages as the drop travels downstream). Along with the prescribed axisymmetry condition on the axis, the boundary conditions are free-slip at the outer radius, free outflow at the channel outlet, and fixed inflow at the channel inlet, where the prescribed incoming gas velocity $U_{0}$ sets the velocity unit. The drop center is initially located on the axis at distance $d=3$ from the inlet boundary, and the velocity field is initially set to zero over the entire domain. Refinement of the adaptive quadtree mesh was allowed up to $400$ cells per drop radius. It is important to appreciate that the choice of an axisymmetric model filters out the three-dimensional instability mechanisms eventually leading to droplet atomisation, and that the observed breakup occurs essentially once the thickness of a liquid sheet becomes smaller than the minimal grid size. Prescribing the same maximum level of refinement for all simulations (except for convergence checks) therefore provides a unique criterion for numerical breakup throughout the study.

The numerical problem can be entirely described using four independent dimensionless parameters: in addition to the Weber number We and the density contrast ${\tt r_{d}}=\rho_{L}/\rho_{G}$ defined in the introduction, we define here the Reynolds number as $\text{Re}=\rho_{G}U_{0}R_{0}/\mu_{G}$ and the viscosity contrast ${\tt r_{v}}=\mu_{L}/\mu_{G}$, where $\mu_{L}$ and $\mu_{G}$ are respectively the dynamical viscosities of the liquid and gas.

I.2 Concerning initial conditions

Dealing with the flow initialisation requires particular care, whether considering our initial velocity field (with $u=0$ over the flow domain and $u=U_{0}$ on the inlet boundary) or the alternative (and equivalent) setups where the velocity field is initialised as $u=U_{0}$ (resp. $u=0$) in the liquid drop and $u=0$ (resp. $u=U_{0}$) in the surrounding gas (with respectively free inflow boundary condition and fixed inlet velocity). It should be emphasized that the numerical model aims at reproducing a physical problem (the impulsive acceleration of a drop traveling across a shock-wave) which as such is a compressible one while discarding the timescale of the acoustic waves. Therefore it necessarily takes a timestep for the numerical solution to adjust and for a (nearly) dipolar velocity field to set in around the drop, which satisfies both the incompressibility constraint and conservation of momentum. As a result however, the relative velocity measured after one timestep between the drop centroid and the gas in the far-field has jumped from the prescribed value $U_{0}$ by an offset $\delta U$ approximately proportional to $\tt r_{d}^{-1}$, so that the deviation is particularly pronounced for small density contrasts.

Importantly, this effect is not a numerical artefact but rather arises from the conservation equations. In the reference frame moving with the shock-wave the drop is initialised with velocity $U_{*}$ and the gas is at rest. As the shock-wave sweeps over the domain a dipolar flow settles in around the spherical drop, which in turn releases a fraction of its momentum as it accelerates the surrounding gas. The overall momentum is conserved: writing ${\cal{V}}$ the volume of the drop, and $U_{1}$ the new relative velocity between the drop and the ambient gas after one timestep, we obtain

[TABLE]

where $C_{M}=\tfrac{1}{2}$ is the added-mass coefficient for a spherical drop magnaudet00 , and we assume the volume of gas entrained by the drop to be equivalent to ${\cal{V}}$. The new relative velocity $U_{1}$ is then

[TABLE]

so that for $r_{d}\gg 1$ the velocity offset $\delta U=U_{1}-U_{*}$ becomes at leading order $\delta U\sim-\tfrac{1}{2}{\tt r_{d}}^{-1}{U_{1}}$.

If we assume that the meaningful quantity is the effective velocity difference between the drop centroid and the far-field resulting from the interaction with the shock-wave, the flow parameters used in previous numerical studies of secondary atomisation might be reconsidered with some caution. Specifically to overcome this issue, our simulations were performed using corrected control parameters ${\text{Re}}_{*},\text{We}_{*}$ built on the anticipated value $U_{*}$, which we determine by requiring the effective velocity difference $U_{1}$ be kept fixed for all the density contrasts considered, and equal to $U_{0}$. The timeseries for the centroid velocity shown in section II are rescaled accordingly, and all the results presented here therefore correspond to the same, consistent definition of the effective We and Re numbers.

II Results

II.1 Deformation at high density contrasts ; comparison with bag breakup experiments

The early stages of drop deformation at large $\tt r_{d}$ are characterised by radial stretching and flattening in the streamwise direction. The drop fate (whether fragmentation or not) then depends on the flow conditions: the restoring effect of surface tension can either overcome the stretching - in which case the drop undergoes oscillatory deformation around a spherical shape while traveling downstream-, or the drop stretches further, inflates and eventually breaks up, displaying a diversity of fragmentation regimes. In the two last decades a vast literature has been dedicated to the experimental investigation of these various deformation modes. As argued in the introduction, most of these experiments were performed in conditions such that the density contrast was typically larger than $500$ and the viscosity contrast $\tt r_{v}$ typically in the $100-1000$ range. In this context, the selection of a particular deformation mode was found to be determined by the Weber number Hinze55 ; ranger69 ; Gelfand74 ; krzeczkowski80 ; wierzba90 ; cao07 ; zhao2010 , with additional influence of the Ohnesorge number whenever the latter exceeds 0.1. At low Ohnesorge numbers (i.e. when viscous effects are weak at the drop scale compared to capillary ones), these studies have shown that the drop dynamics transition from a vibrational mode at the lowest We numbers, where the accelerated drop undergoes capillary oscillations with or without breakup, to so-called bag breakup at $6\lesssim\text{We}\lesssim 15$, where the drop inflates and bursts like a blown balloon, and shear breakup at higher We, where strings of droplets are emitted from the thinning drop edge pilch87 .

Additional distinctions have been introduced between ‘simple’ bag, ‘multiple’ bag or ‘bag-and-stamen’ modes pilch87 , but also between ‘backward-facing’ or ‘forward-facing’ bag modes, some of these features being observable only transiently, and therefore difficult to clearly discriminate from each other. Furthermore, the distinction between ‘simple’ or ‘multiple’ bag deformation can be biased in numerical studies by the neglect of 3D flow features as shown by jain15 . In point of fact such a regime diversity may seem unnecessary from a comprehensive point of view, and for the purpose of the present paper we will restrict ourselves to the distinction between three generic regimes, referred to as oscillatory, bursting or stripping, which are expected to eventually determine the child droplets size distribution. Specifically, the bursting regime is characterised by the formation of a well-defined toroidal rim at the edge of the flattening drop and the inflation and rupture of the central liquid sheet, whatever the curvature of the sheet and the location of its first rupture point (here bursting regime therefore includes any type of bag breakup, whether ‘simple’, ‘multiple’, ‘backward-’ or ‘forward-facing bag’). The stripping regime on the other hand corresponds to the progressive mass loss associated with small droplets being emitted away from the thinning edge (this regime has been also referred to as ‘sheet thinning’ or ‘shear breakup’ in the literature). Finally, oscillatory regime here denotes the absence of any fragmentation process. Note that this attempted classification builds on a phenomenological description of the deformation and breakup dynamics and by no means on an assumed distinction between the physical mechanisms underlying fragmentation. The mechanisms for sheet piercing, whether associated with aerodynamic effects such as the development of a flapping instability, an acceleration-driven hydrodynamic instability, or even the presence of impurities in the drop liquid, are poorly understood and remain an important subject of investigation.

Figure 1 illustrates the drop typical morphological evolution corresponding to these three regimes, as can be observed successively by increasing We at fixed {${\tt r_{d}}=1110$, ${\tt r_{v}}=90.9$, $\text{Re}=1090$} in our simulations, where the contrasts between fluid properties ${\tt r_{d}}$ and ${\tt r_{v}}$ are typically close to that of an air-water system at $293~{}{\rm K}$. The two examples of oscillatory deformation shown in figure 1(a) and (b) respectively correspond to the $\text{We}=2.5$ and $\text{We}=5.65$ cases: the drop flattens and a toroidal rim forms at the edge. Here the case $\text{We}=5.65$ (b) transiently appears to be on the verge of break-up until surface tension effects induce the rim to drain back into the drop core region. The two examples of bursting behaviour in figure 1(c) and (d) respectively correspond to what is classically referred to as bag (or backward-facing bag) mode and multiple bag (or even bag-and-stamen) mode, here for $\text{We}=7.5$ and $\text{We}=12$. Both evolutions are characterised by the formation of a well-defined rim at the edge of the flattening drop and the extreme thinning of the liquid sheet in the interior, eventually causing the inflated drop to burst. Bursting occurs as the bag thickness locally becomes smaller than the minimal grid size, and is typically initiated in the vicinity of the rim neck. It has been argued that the selection of a simple or multiple bag deformation mode could be associated to different most unstable wavelengths with respect to Rayleigh-Taylor instability theofanous04 . In point of fact, whether ‘simple’ or ‘multiple’ bag is seen here to result from an inertia competition between the outer, heavy region forming the rim and the inner sheet at the drop center, resulting in one region being accelerated faster than the other. The size and retractation speed of the rim Taylor59c ; Culick60 being controlled by the surface tension, it is naturally expected that larger We (hence smaller, lighter rims) will result in the rim neck moving outward and the inner region becoming heavier than the rim.

Finally, two stripping-mode examples are shown in figure 1(e) and (f) respectively for $\text{We}=20$ and $\text{We}=40$: here the rim is hardly (if at all) visible, and the early stages of drop deformation exhibit a typical ‘axe-like’ shape (as in the snapshots shown at $\tilde{t}~\sim 40$). This particular feature progressively becomes more acute for larger We as the smoothing operated by the weaker surface tension cannot counteract the erosion of the drop by its own vortical wake (as illustrated by the vorticity fields in figure 1 for $\{\text{We}=40\,;\,\text{Re}=1100\,;\,{\tt r_{d}}=1110\,;\,{\tt r_{v}}=90.9\}$ ; bottom line). As the drop further flattens the indentation carved by the vortical wake tightens and shuts down, until its remnant becomes the noticeable irregularity of the drop edge shown at $\tilde{t}=60$ in figure 1(f). Fragmentation eventually occurs as small droplets are stripped away repeatedly from the thinning edge of the dramatically stretching drop, which behaves like a downstream-facing bag. The transition from bursting to stripping mode appears to be a continuous one and corresponds to the disappearance of the rim.

Figure 2(left) shows the displacement of the drop tip along the streamwise direction (thick solid line) and the drop radius (thin solid line) as a function of the dimensionless time $\tilde{t}=t\,U_{0}/R_{0}$. Here the particular fluid properties ($\tt r_{d},r_{v}$) and flow conditions ($\text{Re},\text{We}$) are chosen so as to reproduce the experimental conditions reported in opfer14 (see their Figure 4, case e-3). Although some differences arise from the comparison between their experimental results and our numerics (for instance the drop tip position slightly displaces upstream during the first flattening stage, and the slope of the drop radius evolution is close to zero in the numerical results at $\tilde{t}=0^{+}$; both features are missing in the experimental curves), these discrepancies can be traced back to the difference in flow initialisation. Indeed, in the experiments the drop is introduced sideways into the channel and has to travel across the boundary layer where it already feels (differential) acceleration by the gas stream before reaching the region of uniform, maximal flow intensity. We therefore expect the transient dynamics to differ from what is monitored in our simulations. Nevertheless, and despite the uncertainties related to initialisation, we find the evolution of the drop radius (dotted line) to agree well with experimental measurements (denoted by the blue stars ; the method used to superimpose the experimental data from opfer14 with our numerics is carefully detailed in the caption of figure 2). Both timeseries display a clear transition between a first deformation stage, in which the maximal drop radius grows (almost) linearly between $\tilde{t}\sim 50$ and $\tilde{t}\sim 110$, and a later stage of exponential growth, although the latter is initiated slightly earlier in the experiments ($\tilde{t}\sim 120$) than in the numerics ($\tilde{t}\sim 130$, corresponding to snapshot ‘e’ in 2). The snapshots (‘f’ to ‘h’) show that this later stage corresponds to the catastrophic inflation of the bag, and the maximal radius is no longer that of the outer rim. (Closer observation of the drop radius timeseries and its successive inflexion points during the first deformation stage (up to $\tilde{t}\sim 110$) actually points towards a behaviour similar to that of a half-oscillation (with dimensionless period $\sim 200$), interrupted at later times by a catastrophic inflation stage. This point will be further discussed later on.)

Numerical and experimental measurements of the tip displacement are more difficult to compare. Although the trends are qualitatively similar, the acceleration of the drop tip seems to be achieved earlier in the experiments - but here again it should be emphasized that in this case the drop has been already progressively accelerated during its travel toward the center of the channel, making the comparison uncertain. An important outcome of our numerics is that they provide for the first time a direct numerical estimate of the drop thickness (solid line in figure 2(bottom, right)), here measured on the symmetry axis. From $\tilde{t}\sim 60$ (which approximately corresponds to the formation of the rim (snapshot ‘c’) up to $\tilde{t}\sim 110$ (snapshot ‘d’), the time evolution of the thickness on the axis (which for this time window is also the minimal thickness) is well represented by an exponential decay (dashed line). The film thickening observed at later times corresponds to a stage where the minimal film thickness is no longer achieved on the axis. Finally, the later stage of drop deformation is characterised by an exponential growth of the bag length $L$ (figure 2(bottom, right, dotted-dashed line), starting around $\tilde{t}\sim 110$ (snapshot ‘d’). As seen from figure 2(left), this time coincides with the end of the first stage in the evolution of the drop radius, as catastrophic bag inflation takes over the ‘half-oscillatory’ behaviour. Breakup occurs (numerically) as the (axisymmetric) sheet minimal thickness reaches the smallest grid size: assuming an initial drop radius of $R_{0}=1.35\text{mm}$ as in opfer14 , this critical thickness would correspond here to $3.2{\mu}m$ (note that in opfer14 bursting is observed experimentally for a thickness of $\sim 1{\mu}m$).

II.2 From large to small density contrasts ; from bursting to stripping

Importantly, the evolution of drop deformation for increasing Weber numbers at small density contrast seems to significantly differ from the large density contrast case, as already observed from early drop breakup simulations zaleski95 ; han01 . Figure 3(left) illustrates this evolution for ${\tt r_{d}}=10$ and We in the $7.5-40$ range, while the Reynolds number and viscosity contrast are chosen as in figure 1. These results are in excellent agreement with previous numerical results of zaleski95 ; han01 ; gaurav2018 . Three observations are readily made. Firstly, we find the transition from oscillatory to fragmentation regime to occur at significantly higher We compared to the large-density series in figure 1. Secondly, for all the We numbers presented here the drop rapidly bends in the downstream direction and develops (at least transiently) into a downstream-facing bag, which - unlike that of figure 1(e) and (f) at high density contrast - exhibits a well-defined rim. No stripping regime could be observed up to $\text{We}=40$, where the inflated drop deforms into a typical jellyfish shape. Thirdly, although a large (and obviously heavy) rim is observed also at large We, the thin (and obviously lighter) inner region seems to lag behind the rim until it either bursts (at high We) or catches up under effect of surface tension at lower We, so that the competition of inertia between the rim and the inner region discussed for the large density contrast case may not seem to apply here. In fact, this apparent contradiction is removed by considering the instantaneous centroid velocity $U_{C}$ of the drop on figure 3(right), where $U_{C}=1/{\cal V}\int_{\cal V}ud\cal V$. (Note that $U_{C}$ is expected to differ from the velocity of the rim or that of the drop apex but is considered here for simplification.) By the time $\tilde{t}=10$ the centroid velocity has already temporary reached a (near to) saturation plateau and the drop does not feel the acceleration at this stage. (The mechanism for temporary stabilisation of the centroid velocity lies in the nonlinear coupling between drop stretching and acceleration, as will be transparent from the simple model developed and discussed later on in section II.3.)

Finally, the transition between deformations regimes observed at low and high density contrasts (respectively in early numerical studies and experimental ones ; see for example the references given in Introduction) is illustrated by the snapshots series in figure 4, where the Weber number was kept fixed at a small value $\text{We}=7.5$ and the density contrast varied for the first time in the $10-2000$ range: the downstream bending of the flattened drop at early times progressively disappears as the density contrast increases, and is no longer observable from ${\tt r_{d}}\sim 100-250$. The times corresponding to the various snapshots in figure 4 are expressed in terms of the reduced (dimensionless) time $t^{*}={\tt r_{d}}^{-1/2}\tilde{t}$.

A particularly remarkable feature in the drop deformation dynamics is the reduced time $t^{*}_{\times}$ at which the rim starts to develop, which in first approximation appears to be almost invariant with the density contrast (here $t^{*}_{\times}\sim 1.8$ was found for $\text{We}=7.5$). This observation can be simply explained by considering the force balance: discarding viscous effects, a rim forms at the edge of the stretching drop when the restoring effect of surface tension overcomes that of inertia, i.e. when

[TABLE]

where $\kappa$ is the curvature of the rim radius as defined on the schematic drawing in figure 5. At this stage, the rim radius is equivalent to the drop half-thickness $e/2$, so that the previous condition rewrites as

[TABLE]

meaning that the Weber number built on the stretching velocity and the drop half-thickness becomes close to unity. Continuity of stagnation pressure across the interface implies that the stretching velocity is related at short times to the gas velocity through $\dot{R}=U_{0}\sqrt{\rho_{G}}/(\sqrt{\rho_{L}}+\sqrt{\rho_{G}})$ - referred to as Dimotakis velocity dimotakis86 -, reducing to $\dot{R}=U_{0}{\tt r_{d}}^{-\tfrac{1}{2}}$ for $\rho_{G}\ll\rho_{L}$. Using mass conservation in the flattening drop, condition (4) thus becomes

[TABLE]

and finally

[TABLE]

yielding $t^{*}_{\times}\sim 1.2$ for $\text{We}=7.5$. Note that this estimate builds on the assumption that $\rho_{G}\ll\rho_{L}$ and therefore fails as $\rho_{G}$ and $\rho_{L}$ become comparable.

The transition from bursting to stripping regime corresponds to the progressive disappearance of the rim. It is associated with the displacement of the most-thinning region from the interior of the sheet (in the vicinity of the rim neck in bursting regime) toward the outer edge (in stripping regime). Since the argument above provides a typical timescale for the formation of the rim, it is therefore natural to compare $t^{*}_{\times}$ with the typical timescale for viscous entrainment in order to predict the transition between the two regimes. As in kekesi14 , the dimensional ‘shear’ breakup timescale $t_{S}\sim e/U_{L}$ (with $U_{L}$ the velocity close to the interface inside the liquid drop) can be estimated using the continuity of viscous stress across the interface:

[TABLE]

where $\delta\sim\sqrt{\mu_{G}R_{0}/\rho_{G}U_{0}}$ is the typical thickness of the laminar boundary layer in the gas blasius08 , and for simplicity $e$ is assumed to provide a relevant lengthscale for the laminar boundary layer in the liquid. The reduced, dimensionless shear breakup timescale is therefore approximated by the relationship

[TABLE]

In kekesi14 , the argument suggested to predict the transition from bursting to stripping regime is equivalent to comparing the typical reduced (dimensionless) time for shear breakup $t^{*}_{S}$ with unity, or comparing its dimensional counterpart with the typical stretching timescale, built on Dimotakis velocity and the drop undeformed radius. However, this approach was motivated by a numerical study performed at constant We and does not take the surface tension into account. Here we compare instead the typical reduced timescale on which shear stresses operate, $t^{*}_{S}$, with the typical reduced timescale for rim formation $t^{*}_{\times}$: even though the latter is (approximately) independent on the density contrast, the former is not, so that the We number characterising the transition between bursting and stripping regimes is found to scale like

[TABLE]

where $\lambda$ is the prefactor inherited from (8). The predicted transition threshold does indeed decrease with increasing density contrasts, consistently with our observations; moreover, this threshold asymptotically tends toward a constant value at large $\tt r_{d}$, in agreement with the phase diagram presented later on in figure 6. We also note that Hsiang95 suggest a prediction $\text{We}\sim\sqrt{Re}$ for a somewhat related transition from “bowl-” to “dome”-shape regime. This prediction however builds on different theoretical considerations to which the authors do not subscribe; in particular, the authors  of Hsiang95  use the undeformed drop radius $R_{0}$ throughout their analysis, whereas we model the variation of the radius $R(t)$ in time.

II.3 Fragmentation threshold

As argued in the Introduction, determining the effect of the liquid-gas density contrast on the fragmentation threshold at low Ohnesorge number has important implications e.g. for CFD modeling purposes, but has not been investigated so far due to the missing overlap between numerically accessible parameters and explored experimental conditions. Taking advantage of Gerris’s capability of accurately handling large $\tt r_{d}$, the influence of the density contrast on the critical We number for fragmentation was characterised for $\tt r_{d}$ varying by nearly 3 orders of magnitude.

Figure 6 presents the phase diagram obtained for fixed gas Reynolds number $\text{Re}=1090$ and viscosity ratio ${\tt r_{v}}=90.9$ in the ${\tt r_{d}}-\text{We}$ plane, where empty symbols denote oscillatory mode (no breakup) and full symbols denote atomisation. Among the latter, discs correspond to bursting and losanges to stripping regime. Mingled symbols on the other hand correspond to a transitional regime where the drop first undergoes dramatic edge-thinning (thus similar to stripping regime dynamics) before a rim starts to form on the outer edge, and breakup first occurs at a thinner region in the sheet interior. Our most important finding here is that the fragmentation threshold drops by one order of magnitude as the density contrast increases from ${\tt r_{d}}=10$ to ${\tt r_{d}}=250$. In the range ${\tt r_{d}}=250-1000$ the critical value lies between $2.5$ and $5.65$, consistently with the prediction of $\text{We}_{c}\sim 3$ formulated by VB09 in the case of a air-water (rain) system (${\tt r_{d}}\sim 1000$). At larger contrast (${\tt r_{d}}=1000-2000$) the fragmentation threshold is found to rise again, although by a less significant amount (a factor two).

This dependency of the critical We number on the density ratio is not predicted by the existing linear models for the evolution of the drop radius Taylor63 ; VB09 , which for low Ohnesorge numbers predict fragmentation for $\text{We}>{\text{We}}_{c}$ regardless of the density contrast. Further, the exponential behaviour of the drop radius for $\text{We}>{\text{We}}_{c}$ theoretically predicted by such models hardly describes the early-time evolution of the drop maximal radius monitored in our simulations (as can be seen for example on figure 8(right)).

However, a key element to understand the decrease in the fragmentation threshold from small to large density contrast (${\tt r_{d}}=10-250$) lies in the loss of the reduced time $t^{*}$-universality which we observe in the evolution of the (streamwise) centroid velocity $U_{C}$, shown in figure 8(left) for constant $\{\text{We}=7.5;\,\text{Re}=1090$;   ${\tt r_{v}}=90.9\}$, and for ${\tt r_{d}}$ varying in the $10-2000$ range. Whereas $t^{*}$ remains the relevant, typical timescale describing the acceleration of the centroid, figure 8(a) shows a non-ambiguous tendency for the drop centroid velocity to increase faster and reach sooner a (temporary) saturation plateau at smaller density contrasts. Interestingly, numerical results of meng01 , who used a compressible flow model at infinite Re and We to study drop breakup behind a shock-wave, are consistent with this observation despite important differences in flow conditions and modeling. Indeed their results for the evolution of the drop centroid velocity as a function of the reduced time $t^{*}$ show a monotonic increase of the early centroid acceleration for increasingly high Mach numbers (figure 14 in meng01 ), i.e. for increasingly high gas density behind the shockwave (following Rankine-Hugoniot jump conditions) and thus increasingly small density contrasts.

The corresponding evolutions of the drop radius $R$ are shown in figure 8(right), with a blowup on early times in the figure inset: as in figure 8(left), interrupted timeseries for ${\tt r_{d}}>40$ denote (numerical) bursting. For all the density contrasts, and whether fragmentation eventually occurs or not, the evolution of the drop radius is similar up to $t^{*}\sim 1$: except for the two largest contrasts (${\tt r_{d}}=1000-2000$) that yield slightly slower drop stretching, the different timeseries (almost) collapse onto a single curve once expressed in terms of the reduced time $t^{*}$. After this first stretching stage, the successive inflection points displayed by the various timeseries suggest that the drop radius initiates oscillations with growing period and amplitude as the density contrast increases from ${\tt r_{d}}=10$ up to (at least) ${\tt r_{d}}=250$, which in the cases where fragmentation occurs (here for ${\tt r_{d}}>40$) are finally interrupted by a new growth and rapid bursting as in figure 2.

In fact, we can show that the smaller inertia of the lighter drops accounts both for the stabilising effect of small density contrasts and for the (non-exponential) time-evolution of the radius monitored in our simulations. To that effect, we assimilate the flattening drop with a disk of uniform thickness $e$, radius $R$ and constant volume ${\cal V}=4/3\pi R_{0}^{3}$, and use cylindrical coordinates to describe the fluid axisymmetric motion. Conservation of momentum in the streamwise direction can be approximated by

[TABLE]

where $C_{d}$ is the drag coefficient determined by the drop shape, so that the centroid acceleration evolves as the square of the drop radius (and vanishes as soon as the drop has reached the gas stream velocity). Note that $C_{d}$ is in truth a time-dependent quantity and is expected to increase as the drop deforms (from a near to spherical shape to a near to cylindrical one), resulting in a faster increase in the centroid velocity; for simplicity we will however assume that $C_{d}$ is constant in the present model.

Using now cylindrical coordinates to describe the fluid axisymmetric motion, conservation of mass inside the expanding disk yields:

[TABLE]

where $\bar{u}_{r}(r,t)=\int_{-e/2}^{e/2}u_{r}(r,z,t)\,dz$ is the radial velocity flux (per unit length), hence $\bar{u}_{r}=r\,e\dot{R}/R$. Following the approach of VB09 , we now approximate the pressure difference in the gas between the high pressure region at the stagnation point on the drop pole ($r=0$) and the lower pressure region at the drop equator ($r=R$) by that of a hyperbolic stagnation flow with stretching rate $\gamma=\alpha(U_{0}-U_{C})/R_{0}$, where $\alpha$ is a shape coefficient, yielding $p_{G}(R)-p_{G}(0)=-\rho_{G}\gamma^{2}R^{2}/8$. Importantly, we will however assume here that the relative velocity $U_{0}-U_{C}$ remains time-dependent.

The pressure difference in the liquid is estimated from the pressure difference in the gas: at the pole $r=0$, the pressure jump across the interface vanishes due to the negligible curvature of the latter. At the equator however, this pressure jump can be approximated by $2\sigma/e\sim\tfrac{3}{2}\sigma R^{2}/R_{0}^{3}$ according to Laplace’s law, so that the pressure difference in the liquid across the drop radius is finally

[TABLE]

Discarding viscous damping, the axisymmetric Euler equation for the radial velocity $u_{r}$ inside the drop writes

[TABLE]

Assuming uniform radial motion through the drop thickness (hence $u_{r}=\bar{u}_{r}/e=r\dot{R}/R$) and averaging (13) over the drop domain using (12) yields

[TABLE]

an evolution equation for the drop radius non-linearly coupled with that of the centroid velocity derived from (10):

[TABLE]

with $R(0)=R_{0}$ and $U_{C}(0)=0$. The non-linear coupling between (14) and (15) accounts for the role played by the density contrast in the fragmentation threshold. If we ignore the drop acceleration as in VB09 , as is realistic for large $\tt r_{d}$ at early times, (14) becomes linear in $R$. The drop radius then undergoes exponential stretching as soon as the Weber number exceeds the critical value ${\text{We}}_{c}=12/\alpha^{2}$, on a typical timescale determined by the reduced time $t^{*}={\tt r_{d}}^{-1/2}\tilde{t}$. (In VB09 the shape coefficient determining the stretching rate of the hyperbolic stagnation flow in the gas is $\alpha=2$, an intermediate value between the case of a spherical obstacle ($\alpha=3$) and that of a disc ($\alpha=\pi/2$). This value of $\alpha$ yields ${\text{We}}_{c}=3$, in excellent agreement with the fragmentation threshold found numerically at large density contrast.) However, taking (15) into account modifies the critical Weber number: it is clear from (14) that an increase in $U_{C}$ achieved sufficiently fast over the reduced timescale can significantly raise the fragmentation threshold. The acceleration at early times ($\tilde{t}=0^{+}$) is larger for smaller density contrasts, so that a faster decrease in the right-hand side of (14) is expected as $\tt r_{d}$ decreases, associated with a smaller stretching amplitude and oscillation period (expressed in reduced timescale).

Solutions of (14)-(15) determined by numerical integration are shown in figure 8 for constant $\text{We}=7.5$ and three different density contrasts in the $10-1000$ range. Even though our simple non-linear model does not allow for a precise quantitative description of the drop deformation dynamics and does not account for the slight increase in $\text{We}_{c}$ numerically observed at ${\tt r_{d}}=1000-2000$ (presumably because of the important simplifications made in the derivation of (14)-(15)), the principal features of figure 8 - stronger acceleration and earlier (temporary) saturation of the drop velocity for small density contrasts, as well as shorter and weaker oscillations of the drop radius - are well reproduced. Our findings thus demonstrate that the variation in the drop inertia indeed accounts for the dependency of the deformation dynamics on the density contrast over the reduced timescale $t^{*}$, and for the stabilising effect of small density contrasts with respect to fragmentation. However, it is important to emphasize that our model does not address the deformation and bursting dynamics during the later inflation stage. Moreover, our axisymmetric simulations cannot address the inherently three-dimensional mechanisms of sheet piercing in the inflating drop, which would require additional (and computationally much more demanding) studies. Full understanding of the bursting mechanisms, including the possible development of hydrodynamic instabilities in the rim or and the thinning film, is beyond the scope of the present article and remains a challenging problem for future investigation.

III Conclusion

Our axisymmetric simulations of impulsively-accelerated drop breakup show for the first time a significant influence of the liquid-gas density contrast on the fragmentation thresholds at low Ohnesorge number. The density contrast was varied over two decades (${\tt r_{d}}=10-2000$), showing good agreement with previous numerical and experimental works. We suggest a simple theoretical argument to determine the critical Weber number characterising the transition from bursting to stripping regime, which we find to increase for decreasing density contrasts, consistently with our numerical results. Although the objective of the present paper was not to provide a review of the studies that could experimentally assess the validity of these theories over the whole parameter space (large and small Reynolds and Weber numbers, and large and small density/viscosity contrasts), we believe that the full range of such experimental investigations remains to be done and would be of great interest. Finally, we explain the significant rise of the fragmentation threshold observed here at small density contrasts by the nonlinear coupling between the stretching drop radius and the centroid acceleration. Our results demonstrate the need for adaptation of existing breakup criterion, such as those currently used in non-DNS CFD models, so as to take into account the effect of density contrast on the fragmentation threshold and regimes.

Acknowledgement. The authors wish to thank Sasha Korobkin, Stéphane Popinet, Wojciech Aniszewski and Christophe Josserand for fruitful discussions.