J. Fluid Mech. (2018), vol. 835, pp. 1108–1135.
c
Cambridge University Press 2017
doi:10.1017/jfm.2017.804
1108
Numerical simulation of the aerobreakup of a
water droplet
Jomela C. Meng
1,
† and Tim Colonius
1
1
California Institute of Technology, Pasadena, CA 91125, USA
(Received 29 October 2016; revised 28 October 2017; accepted 3 November 2017;
first published online 29 November 2017)
We present a three-dimensional numerical simulation of the aerobreakup of a spherical
water droplet in the flow behind a normal shock wave. The droplet and surrounding
gas flow are simulated using the compressible multicomponent Euler equations
in a finite-volume scheme with shock and interface capturing. The aerobreakup
process is compared with available experimental visualizations. Features of the droplet
deformation and breakup in the stripping breakup regime, as well as descriptions of
the surrounding gas flow, are discussed. Analyses of observed surface instabilities and
a Fourier decomposition of the flow field reveal asymmetrical azimuthal modulations
and broadband instability growth that result in chaotic flow within the wake region.
Key words: breakup/coalescence, drops, shock waves
1. Introduction
The study of droplet aerobreakup has historically been motivated by three
applications: bulk dissemination of liquid agents, raindrop damage during supersonic
flight, and secondary atomization of liquid jets in turbomachinery. Much of the
aerobreakup literature has focused research efforts on characterizing and mapping
various breakup regimes (e.g. Lane 1951; Engel 1958; Hanson, Domich & Adams
1963; Ranger & Nicholls 1968; Hsiang & Faeth 1995), calculating characteristic
breakup times (e.g. Ranger & Nicholls 1968; Hsiang & Faeth 1992), quantifying
dependence on parameters such as density and viscosity ratios (e.g. Hanson et al.
1963; Theofanous et al. 2012), predicting final drop size distributions (e.g. Ranger
& Nicholls 1968; Pilch & Erdman 1987), and quantifying unsteady drag properties
(e.g. Engel 1958; Simpkins & Bales 1972; Joseph, Belanger & Beavers 1999).
Unfortunately, these experimental and theoretical research efforts have resulted in
many, and often conflicting, phenomenological models describing the aerobreakup
process, and to date, a definitive understanding of aerobreakup remains elusive
(Khosla, Smith & Throckmorton 2006).
Beginning with the work of Hinze (1949), the Weber number, We = ρ
g
u
2
g
D
0
/σ , has
been the principal parameter used to delineate the various regimes of aerobreakup.
Traditionally, there exist five distinct regimes that are well established in the literature
(Pilch & Erdman 1987; Guildenbecher, López-Rivera & Sojka 2009). They are,
† Email address for correspondence: jomela.meng@caltech.edu
Downloaded from https://www.cambridge.org/core. Caltech Library, on 20 Dec 2017 at 20:54:57, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2017.804
Numerical simulation of the aerobreakup of a water droplet 1109
in order of increasing We, the vibrational, bag, bag-and-stamen, stripping, and
catastrophic regimes.
The stripping regime, which is of primary interest in this paper, marks a transition
in breakup physics that fundamentally differs from that of the preceding breakup
regimes. Generally speaking, the stripping regime is characterized by an initial
deformation of the droplet into a disk-like shape. Following the shape change, droplet
fluid is observed to be stripped from the droplet’s periphery in a region near the
droplet equator (defined as a polar, or inclination, angle of ϕ = π/2). Ranger &
Nicholls (1968) postulated the ‘boundary layer stripping’ or ‘shear stripping’ model,
where boundary layers both inside and outside the drop become unstable at the
droplet equator, and are subsequently stripped off by the ambient gas flow. Liu &
Reitz (1997) later proposed an alternate mechanism known as ‘sheet thinning’, where
the droplet is initially flattened by the pressure gradient between the drop’s poles
(ϕ = 0, π) and equator. Once flattened, the strong inertial forces from the surrounding
flow draw a thin sheet of liquid off the periphery. This sheet is accelerated, stretched,
and bent in the direction of flow, and eventually breaks up into streamwise ligaments
that fragment into individual drops. Due to the flattened disk-like shape of the
deformed droplet, flow separation occurs for all practical values of the Reynolds
number, Re = ρ
g
u
g
D
0
/µ
g
, and the sheet thinning mechanism can be considered an
inviscid phenomenon with no dependence on Re (Guildenbecher et al. 2009).
The instability of the thin liquid sheet that is drawn from the droplet’s periphery in
the sheet thinning model is thought to be responsible for the generation of product
droplets. Liu & Reitz (1997) described a ‘stretched streamwise ligament breakup’
mechanism wherein, for low liquid flow rates of a planar liquid sheet sandwiched
between two shear air layers, streamwise vortical waves, alternating with thin liquid
membranes, would grow along the sheet. The membranes would burst first due to the
rotation of the streamwise vortices, leaving streamwise ligaments that subsequently
broke up. More recently, Jalaal & Mehravaran (2014) attributed the breakup to the
rise of Rayleigh–Taylor (RT) instability waves on the sheet. In this mechanism, the
Kelvin–Helmholtz (KH) instability generates axisymmetric waves at the liquid–gas
interface. The transient acceleration of these wave crests or rims (in the case of
the liquid jet) into the downstream air triggers a RT instability, which produces
‘transverse azimuthal modulations.’
Recent work by Theofanous, Li & Dinh (2004) studying aerobreakup in rarefied
supersonic flows, and subsequent publications (Theofanous & Li 2008; Theofanous
2011; Theofanous et al. 2012), has substantially changed the overall understanding of
aerobreakup. Theofanous et al. (2004) proposed a reclassification into two principal
breakup regimes: Rayleigh–Taylor piercing (RTP) and shear-induced entrainment (SIE).
Theofanous & Li (2008) described SIE as a combination of shear-driven radial motion,
which results in the flattening, as well as instabilities on the stretched liquid sheet.
Perhaps most importantly, this reclassification argues that the catastrophic breakup
regime does not exist. Theofanous & Li (2008) contended that the wavy interface
on the upstream side of the droplet was an artefact created by a projected view of a
complex flow field, and that no RT waves exist or pierce the drop. Using laser-induced
fluorescence as their experimental visualization technique, SIE was proposed as the
terminal regime for We > 10
3
.
Recently, numerical simulation has emerged as a tool for studying aerobreakup.
Unfortunately, due to the high computational costs of fully three-dimensional (3D)
simulations, numerical aerobreakup studies have often invoked two-dimensional
(2D) (Zaleski, Li & Succi 1995; Igra & Takayama 2001a,b,c; Chen 2008) or
Downloaded from https://www.cambridge.org/core. Caltech Library, on 20 Dec 2017 at 20:54:57, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2017.804
1110 J. C. Meng and T. Colonius
axisymmetric (Han & Tryggvason 2001; Aalburg, Leer & Faeth 2003; Wadhwa, Magi
& Abraham 2007; Chang, Deng & Theofanous 2013) approximations. Additionally,
fluid density ratios are often assumed to be small, or the fluids are considered to
be incompressible (Han & Tryggvason 2001; Aalburg et al. 2003; Quan & Schmidt
2006; Jalaal & Mehravaran 2014; Xiao, Dianat & McGuirk 2014; Castrillon Escobar
et al. 2015; Jain et al. 2015). Previous work from the authors (Meng & Colonius
2015) simulated the early stages of two-dimensional aerobreakup, with comparison
to the experimental work of Igra & Takayama (2001c). Qualitative features of the
breakup process, such as the presence of a transitory equatorial recirculation region
and an upstream jet in the wake, were discussed. Additionally, a parametric study
varying incident shock strength was performed to study the effects of the transition
between subsonic and supersonic post-shock flow. A novel method of calculating
the cylinder’s centre-of-mass properties was also utilized to obtain accurate unsteady
acceleration and drag histories.
Using numerical simulation, the purpose of this paper is to elucidate the physical
breakup mechanisms responsible for the fully three-dimensional, compressible
aerobreakup of a single water droplet in air. While a direct numerical simulation
would be ideal for such an investigation, the computational grid required to fully
resolve all scales of the aerobreakup problem is intractable on currently available
computational resources; a compromise must be made. We thus consider the ‘inviscid’
case, which is in turn regularized by artificial viscosity and numerical diffusion
of the interface. The effects of these approximations are discussed in detail in
§ 3.2. Additionally, we will attempt to interpret the numerical results in a manner
consistent with the uncertainties introduced by the approximations. These results
fill a gap in the current state of droplet aerobreakup knowledge associated with
the underlying fundamental flow physics that dictates the experimentally observed
phenomena. Building upon the computational efforts of Coralic & Colonius (2014),
we utilize a compressible multicomponent flow solver to numerically investigate this
fundamental fluid dynamics problem. The rest of this paper is organized as follows.
In § 2, we describe the problem set-up and introduce the governing equations. The
numerical algorithm is subsequently described in § 3, along with the calculation of
various flow quantities used in the analysis. We describe the evolution of the liquid
droplet, show its centre-of-mass properties, and compare with available experimental
visualizations in §§ 4, 4.1 and 4.2. The flow field in the gas phase is investigated in
§ 4.3. The mechanisms of surface instabilities are discussed in § 5.1, and are followed
by a Fourier decomposition analysis in § 5.2. Finally, concluding remarks are made
in § 6.
2. Physical modelling
2.1. Problem description
In the laboratory, shock tubes are most often used to generate large relative velocities
between the liquid droplet and surrounding gas flow. Normal shock waves, in and
of themselves, have little effect on the droplet. Instead, they serve as a reliable and
repeatable way to create a high-speed flow around the droplet, which is responsible
for the droplet’s subsequent deformation and disintegration (Hanson et al. 1963;
Ranger & Nicholls 1968; Joseph et al. 1999). In figure 1, a schematic is shown of
the initial condition and computational grid. The shock wave of strength M
s
= 1.47
(modelled as a step discontinuity in fluid properties) is travelling in air towards
the water droplet, which is initialized as a spherical interface with diameter D
0
on
the cylindrical grid. The drop and the ambient air downstream of the shock are
Downloaded from https://www.cambridge.org/core. Caltech Library, on 20 Dec 2017 at 20:54:57, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2017.804
Numerical simulation of the aerobreakup of a water droplet 1111
Post-shock
air
Pre-shock
air
NRBC
NRBC
NRBC
NRBC
r
z
FIGURE 1. (Colour online) Schematic of the 3D initial condition and non-uniform
computational domain at 1:10 of the actual resolution. The grid extends radially outward
for 6D
0
, and the axial extents are −7D
0
6 z 6 15.5D
0
.
initialized at rest. The entire 3D flow field is simulated with no enforced symmetries.
Non-reflective boundary conditions (NRBC), applied at all computational domain
boundaries, approximately extend the surrounding air to infinity. The simulation is
performed on the computational domain Ω = [−7D
0
, 15.5D
0
] × [0, 6D
0
] × [0, 2π]
with a spatial resolution of (N
z
, N
r
, N
θ
) = (800, 600, 320) grid cells. The grid is
stretched towards the axial boundaries using a hyperbolic tangent function. The most
refined portion of the grid is located near the initial location of the droplet and
in the region of the near-field wake. In this region, the nominal axial and radial
grid resolution corresponds to 100 cells per original droplet diameter. The azimuthal
resolution is chosen such that the cells near the spherical droplet interface are close
to regular. Previous testing of grid resolution sensitivities in 2D (Meng & Colonius
2015; Meng 2016) suggests the present spatial resolution captures the salient flow
features without being computationally cumbersome. The simulation is performed with
a constant Courant–Friedrichs–Lewy (CFL) number of 0.2. Since the initial condition
is axisymmetric, the gas is initially seeded with small random radial and azimuthal
velocity perturbations, the largest of which have approximate magnitudes of O(10
−4
u
s
),
where u
s
is the post-shock gas velocity. These white-noise-type perturbations are
generated via the Fortran compiler’s intrinsic random number generator, and are
applied to both the pre-shock and post-shock gas in the initial condition.
2.2. Governing equations
We model the flow with the compressible multicomponent Euler equations. In addition
to being compressible, each fluid is considered immiscible and does not undergo phase
change. In the absence of mass transfer and surface tension, material interfaces are
Downloaded from https://www.cambridge.org/core. Caltech Library, on 20 Dec 2017 at 20:54:57, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2017.804
1112 J. C. Meng and T. Colonius
simply advected by the local flow velocity. For a system of two fluids, this gives
rise to a five-equation model, introduced in its inviscid form by Allaire, Clerc &
Kokh (2002). Implications of the choice to neglect viscous and capillary effects are
addressed in detail in § 3.2. The five-equation model, (2.1), consists of individual
continuity equations for each of the fluids, (2.1a) and (2.1b), mixture momentum,
(2.1c), and energy, (2.1d), equations, and a transport equation for the gas volume
fraction, (2.1e):
∂(α
g
ρ
g
)
∂t
+ ∇ · (α
g
ρ
g
u) = 0, (2.1a)
∂(α
l
ρ
l
)
∂t
+ ∇ · (α
l
ρ
l
u) = 0, (2.1b)
∂(ρu)
∂t
+ ∇ · (ρu ⊗ u + pI) = 0, (2.1c)
∂E
∂t
+ ∇ · ((E + p)u) = 0, (2.1d)
∂α
g
∂t
+ u · ∇α
g
= 0, (2.1e)
where ρ is the density, α is the volume fraction, u is the velocity vector, p is the
pressure, E is the total energy defined as E = ρε + (1/2)ρkuk
2
, ε is the specific
internal energy, and the subscripted g, l represent, respectively, the gas and liquid
fluids. Each equation in the five-equation model, (2.1), is evolved independently of
all others. While alternate models for multicomponent flows exist within the literature
(e.g. Kapila et al. 2001; Murrone & Guillard 2005; Saurel, Petitpas & Berry 2009;
Pelanti & Shyue 2014), the present model is chosen for its simplicity and desirable
conservation properties. The simplified interface model theoretically ensures volume
fraction positivity, and does not require regularization to fully specify shock jump
conditions. Though it is not capable of physically handling mixture regions, the model
is sufficient for the immiscible fluids assumed in the multicomponent flows of interest.
The simulation of material interfaces is made possible by the volume-of-fluid
method, which belongs to the broader class of interface-capturing schemes. One
common characteristic of these schemes is the relaxation of the natural sharpness
of material discontinuities. That is, the interfaces are allowed to numerically diffuse,
resulting in an interface region of small, but finite, thickness. Within this interface
region, a non-physical fluid mixture exists, which is appropriately treated using
mixture rules that are defined in § 2.2.1. For numerical stability purposes, material
interfaces are not initialized as sharp discontinuities, but are smeared over a few grid
cells. Based on previous results (Johnsen 2007; Coralic 2015) and sensitivity testing
by the authors for the specific problem of aerobreakup, this initial artificial smearing
over a few cells is known to have negligible impact on the computational results.
2.2.1. Equation of state
The five-equation model, (2.1), is closed with the specification of an appropriate
equation of state (EOS) that relates the fluid densities, pressures, and internal energies.
The stiffened gas EOS (Harlow & Amsden 1971), is used in our flow solver to model
both gases and liquids:
p = (γ − 1)ρε − γ π
∞
, (2.2)
Downloaded from https://www.cambridge.org/core. Caltech Library, on 20 Dec 2017 at 20:54:57, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2017.804
