by user

Category: Documents





Contact us
My IOPscience
Numerical simulation of droplets, bubbles and waves: state of the art
This article has been downloaded from IOPscience. Please scroll down to see the full text article.
2009 Fluid Dyn. Res. 41 065001
View the table of contents for this issue, or go to the journal homepage for more
Download details:
IP Address:
The article was downloaded on 17/01/2011 at 20:49
Please note that terms and conditions apply.
Fluid Dyn. Res. 41 (2009) 065001 (24pp)
Numerical simulation of droplets, bubbles and
waves: state of the art
Daniel Fuster1,2 , Gilou Agbaglah1,2 , Christophe Josserand1,2 ,
Stéphane Popinet3 and Stéphane Zaleski1,2,4
UMR 7190, Institut Jean Le Rond d’Alembert, Univ. Paris 06, F-75005 Paris, France
CNRS, UMR 7190, Institut Jean Le Rond d’Alembert F-75005 Paris, France
3 National Institute of Water and Atmospheric Research, Wellington, New Zealand
E-mail: [email protected]
Received 5 August 2009, in final form 16 November 2009
Published 4 December 2009
Online at stacks.iop.org/FDR/41/065001
Communicated by M Funakoshi
This work presents current advances in the numerical simulation of twophase flows using a volume of fluid (VOF) method, balanced-force surface
tension and quad/octree adaptive mesh refinement (AMR). The simulations
of the atomization of a liquid sheet, the capillary retraction of a liquid sheet
and two- and three-dimensional wave breaking all for air/water systems are
used to show the potential of the numerical techniques. New simulations of
atomization processes for air/water conditions are allowing us to investigate
the processes leading to the appearance of instabilities in the primary
atomization zone in real conditions. For the retracting liquid sheet, the new
simulations show that two different regimes can be encountered as a function
of the Ohnesorge number. For large values, a laminar flow is encountered
inside the rim and a steady state is reached after a quick transient state. For
small values, a turbulent flow is generated inside the rim, which is responsible
for large oscillations in the rim size and neck thickness. The breaking wave
case study demonstrates the orders-of-magnitude efficiency gains of the AMR
(Some figures in this article are in colour only in the electronic version)
Author to whom any correspondence should be addressed.
© 2009 The Japan Society of Fluid Mechanics and IOP Publishing Ltd
Printed in the UK
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
1. Introduction
The physics and particularly the dynamics of two-phase flow have grown into a major
scientific domain with crucial fundamental issues as well as many practical and industrial
applications. Fields of interest include drop impact phenomena involved in the study of
rain on soils or foliage, inkjet printing and combustion chambers (Yarin and Weiss 1995,
Josserand and Zaleski 2003). In particular, atomization processes (Eggers and Villermaux
2008) play an important role in combustion studies. Two-phase flow is also prevalent in
geophysical contexts such as river motion and meandering or ocean dynamics (Coantic 1980).
The dynamics of interfaces is rich and complex, with many mathematical and technical issues
that need to be resolved to account for the specific role of these interfaces (a surface in three
dimensions (3D), a line in 2D) separating two domains with different physical properties.
Recently, the experimental access to high-speed imaging and the systematic use of self-similar
approaches and matched asymptotics methods (see, for instance, Eggers 1997) have greatly
improved our general understanding of multiphase flows. These developments and the recent
improvements in numerical techniques have also made it possible to tackle problems that were
still inaccessible a few years ago, such as 3D impacts or spray dynamics. All these aspects
have tremendously increased interest in the ever more complex and realistic two-phase flow
simulations that are now feasible. The aim of this paper is thus to present recent achievements
in two-phase flow study based on up-to-date numerical methods. The numerical techniques
used here combine two classical methods, namely the volume of fluid (VOF) method that
allows a good treatment of the interface, and the adaptive mesh refinement (AMR), which
optimizes the numerical calculations as discussed earlier by Popinet (2009) and Fuster et al
(2009). In this paper, after a brief review of important recent developments in numerical
methods, we illustrate what is now achievable through three typical two-phase flow studies in
air–water conditions.
2. Numerical method
Most multiphase flow problems for practical applications exhibit all or several of the following
characteristics: high surface tension, low viscosity, high density ratios, complex and evolving
interface topologies and spatial scales ranging over several orders of magnitude. In this
context, an ideal numerical method for the solution of the two-phase Navier–Stokes equations
with surface tension would have the following properties:
• robust representation of evolving, topologically complex interfaces,
• accurate representation of surface tension, which requires accurate estimates of interface
normal and curvature,
• robust and accurate handling of large density and viscosity ratios,
• efficient representation of evolving flow features of widely different characteristic spatial
While distinct numerical methods have one or several of these properties, arguably
no method published to date possesses all of them together. Methods using an implicit
representation of the interface such as VOF (Scardovelli and Zaleski 1999) or level set (LS)
(Sussman et al 1994) can robustly and efficiently represent evolving, topologically complex
interfaces but generally suffer from an inaccurate representation of surface tension (Popinet
and Zaleski 1999, François et al 2006). Methods using an explicit representation of the
interface such as arbitrary Lagrangian–Eulerian (ALE) discretizations (Fyfe 1988, Dai and
Schmidt 2005) or front-tracking (FT) (Popinet and Zaleski 1999, Shin et al 2005) can provide
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
an accurate representation of surface tension but have difficulty dealing with complex and
evolving interface topologies. Finally, dynamic mesh-adaptive methods can deal efficiently
with phenomena involving a wide range of spatial scales but their implementation, to date,
has been limited to a small subset of the methods referred to above.
Recently, one of us proposed new extensions of these methods that allow us to
combine their respective strengths and thus come closer to the ideal numerical method
described above (Popinet 2009). The basis for the numerical scheme is the incompressible
adaptive Navier–Stokes solver Gerris (Popinet 2003, 2007), which uses an octree for spatial
discretization. This type of discretization was shown to lead to efficient adaptivity for singlephase flows (the cost of adapting the mesh at every time step is generally less than 5% of the
total computing effort). This approach is different from the more classical block-structured
adaptive schemes of Almgren et al (1998) in that each discretization element can be refined
independently from the others. This fine-scale adaptivity was shown in Popinet (2009) to
be very relevant to flows with interfaces where it is often desirable to refine on and along
interfaces only. While doing this is possible with block-AMR techniques, the finite size of
AMR blocks leads to significant increases in total mesh sizes.
A challenge is then to extend classical interface-tracking techniques to these adaptive
meshes. A simple approach is to restrict adaptivity only to the bulk of each phase. This
guarantees that the interface between phases never changes resolution and leads to a trivial
implementation of existing methods (Greaves 2004). This comes at the cost of never being
able to refine along the interface. As was shown by Popinet (2009) and as is further illustrated
in the present paper, this restriction can lead to a dramatic decrease of the efficiency of
adaptivity for realistic problems. To lift this restriction, Popinet (2009) proposed extending
the VOF method to adaptive octree grids. The technique rests on the creation of virtual regular
Cartesian stencils, reconstructed dynamically from the underlying octree discretization. Using
the algorithmic properties of octrees, simple and efficient algorithms can be designed to
achieve this. The resulting technique allows full adaptivity along the interface while keeping
all the other advantages of the VOF method: arbitrary interface topology and topological
changes, good mass conservation properties and sharp interface description. Figures 1 and 17
illustrate VOF-reconstructed interfaces with a spatial resolution varying along the interface
obtained with the Gerris solver.
An oft-mentioned drawback of VOF methods is the difficulty of obtaining accurate values
for useful geometrical properties such as interface normal directions or curvature. This is
often contrasted with the relative simplicity with which these values can be obtained with
e.g. LS techniques. A recent study by François et al (2006) casts doubt on this simplistic
interpretation however and shows that obtaining accurate second-order curvature estimates
in the LS method is not trivial. This paper also compares LS-derived curvature estimates
with a ‘rediscovered’ height-function (HF) curvature estimation using a VOF fraction field.
It shows that the HF scheme gives a second-order accurate estimate of curvature from VOF
volume-fraction fields while being simple to implement on regular Cartesian grids. While this
is true for well-resolved interfaces, a problem occurs when the interface curvature becomes
comparable to the mesh size. In this limit the standard HF method becomes inconsistent.
The problem is addressed in Gerris by switching to a novel generalized HF method for low
resolutions, which ensures consistent second-order convergence of the curvature estimates
across all resolutions (Popinet 2009). Gerris further generalizes the HF method to octree grids
by using the virtual regular Cartesian stencils method mentioned previously.
Although computing accurate curvature estimates is an essential step to obtain an accurate
surface tension formulation, it is not sufficient. The inaccuracy of surface tension formulations
has been the traditional weak point of methods using an implicit representation of the interface
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Figure 1. Jet created by the collapse under surface tension of an air bubble just below the surface
of still water. The interface is represented using the VOF-reconstructed fragments in each cell.
The resolution varies along the interface according to local interface curvature. The finest scales
(droplets and jet) are 128 times more refined than the coarsest scales (flat surface in the corners
of the picture). This 3D, two-phase, air–water simulation can be compared with the axisymmetric,
free-surface (i.e. single phase) simulations of Duchemin et al (2002). Accurate surface tension,
robustness for large density ratios and changes of topology as well as adaptivity are all needed to
obtain these results.
(such as VOF or LS). A classical simple test case is that of a bubble or droplet at rest in a
surrounding fluid. In this case the exact solution is given by Laplace’s law stating that—for a
constant surface tension coefficient—the pressure jump across the interface exactly balances
the surface tension force, resulting in a stationary, spherical interface (Popinet and Zaleski
1999). Until recently, all surface tension formulations coupled with an implicit interface
representation failed this simple test case and produced non-equilibrium solutions with the
so-called ‘spurious’ or ‘parasitic’ currents of varying amplitude depending on parameters
such as viscosity, density ratio, spatial and temporal resolution. At worst, these parasitic
currents can be strong enough to dominate the solution. This is of particular concern for
practical applications, which often fall in the parameter range where parasitic currents are
strong. Although there has been a vast number of papers trying to address the problem in
the past ten years, the method implemented in Gerris and published by Popinet (2009) is
the first that demonstrates recovery of an exact equilibrium solution for this simple test case
(see table 1). Furthermore, this exact solution was also shown to be obtained irrespective of
spatial resolution or physical parameters. The method stresses the importance of the concept
of ‘balanced-force’ continuum formulation introduced by Renardy (2002) and couples it with
the generalized HF method to obtain exact equilibrium solutions.
The overall scheme is simple to implement, only requires a single VOF interface
representation (in contrast to coupled methods such as VOF/markers and VOF/LS), works
on octree grids allowing adaptivity along the interface, conserves mass and provides accurate
surface tension solutions for a range of test cases.
Table 1 summarizes the progress made in tackling the problem of ‘parasitic
currents’ around a stationary spherical droplet with surface tension. The initial VOF
VOF-CSS (Gueyffier et al 1998)
FT-Markers (Popinet and Zaleski 1999)
GFM (Kang et al 2000)
VOF-PROST (Renardy and Renardy 2002)
FT-Markers (Shin et al 2005)
LS-VOF (Sussman et al 2007)
DG (Compére et al 2008)
LS (Herrmann 2008)
VOF-HF (Popinet 2009)
12 000
12 000
728 000
12 000
12 000
12 000
12 000
612 000
Spatial resolution
Table 1. Amplitude of spurious currents for a selection of methods √
published in the literature. The spatial resolution is given in number of grid points per droplet diameter. The amplitude
of spurious currents is given in non-dimensional form as: max |u| D/σ with D the droplet diameter and σ the surface tension coefficient. La is the Laplace number: σρ D/µ2 with
µ being the viscosity and ρ the density of the liquid.
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
continuum–surface–stress implementation of Gueyffier et al (1998) had very strong spurious
currents independently of spatial resolution. This problem was addressed by Popinet
and Zaleski (1999) by switching to an FT (markers connected by third-order splines)
representation of the interface coupled with a ‘surface-tension corrected’ pressure gradient.
While spurious currents were drastically reduced using this approach and showed convergence
with resolution, the marker technique was only 2D and could not deal with topology changes.
The PROST method of Renardy and Renardy (2002) was the first VOF method with low
spurious currents that we are aware of. The important concept of ‘balanced-force’ surface
tension representation in the context of continuum-surface-force methods is introduced. The
paraboloid-fitting method used to compute curvature from volume fractions was highly
computationally expensive, however. Most of the other methods published since then were
refinements of these earlier methods and showed similar amplitudes of spurious currents. The
VOF, height-function, balanced-force method implemented in Gerris is the first technique that
demonstrated convergence to exact equilibrium (to within machine precision) irrespective of
Ohnesorge number and spatial resolution (Popinet 2009).
Finally, large density ratios have also proven difficult to handle using methods with an
implicit representation of the interface. For such methods it is tempting and often necessary
to smooth the jumps in physical properties (including density) occurring at the interface.
While this leads to simple finite-difference discretizations of the (now continuous) partial
differential equations describing two-phase flows, the resulting numerical approximations
can lead to physically inconsistent results. Several schemes have been proposed to try to
address some of these problems: for example, both the ‘pressure gradient correction’ scheme
(Popinet and Zaleski 1999) and the ‘ghost-fluid method’ (Kang et al 2000) introduce better
(and similar) discretizations of the pressure jump due to surface tension forces but do not
address the discontinuity in fluid momentum and/or viscous stresses due to jumps in density
and viscosity. While mass conservation has always been a key advantage of VOF methods
(compared to e.g. LS methods), momentum conservation has generally received less attention.
Many implementations thus discretize the non-conservative form of the variable-density
Navier–Stokes equations, which has the advantage of using the continuous velocity field as a
primary variable, rather than the discontinuous momentum field. The density field and velocity
field are then advected using different schemes (e.g. VOF for density and centered/upwind
scheme for velocity), which can lead to anomalous numerical momentum diffusion across
the interface. To address this problem Rudman (1998) proposed to advect momentum using
the same VOF scheme as for density so that momentum conservation was recovered. This
scheme relies on a sharp interface representation and is not compatible with smoothing. We
intend to implement this scheme in Gerris but at present we use a non-momentum-conserving
velocity/density formulation.
A related problem for large density ratios is that the Poisson problem for the pressure—
derived from the incompressibility condition—becomes stiffer as the density jump increases.
This can again be somewhat minimized by smoothing the interface, but as discussed
previously this is not a satisfactory solution (although this is often done in practice). Gerris
uses the multigrid Poisson solver presented by Popinet (2003) for single-phase flows that fit
nicely with the octree discretization. This solver generally works well also for two-phase
flows with large density ratios (as demonstrated in this paper); however, its performance
(in terms of convergence speed) depends—as for any geometric multigrid solver—on the
consistency of the representation of domain (i.e. interface) topology on successively coarser
grids. For complex flows this consistency is not guaranteed and this can lead to important
performance penalties. A range of recent linear solver techniques could be used to address this
problem, including: algebraic multigrid, multigrid-accelerated GMRES (generalized minimal
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Torres and Brackill (2000)
Herrmann (2008)
Second-order convergence
Frequency error (%)
Number of grid points per droplet diameter
Figure 2. Convergence of the error between the computed and theoretical frequencies for the
mode 2 oscillations of a liquid droplet.
residual method), Krylov-subspace methods, etc and we plan to explore this issue further in
the future.
To study the accuracy of the surface-tension representation in the case of large density
ratios, we considered the case of 2D oscillations of a liquid drop floating in a lighter fluid. The
density ratio is 1/1000. No explicit viscosity is added so that any damping of the oscillations
is due to numerical diffusion only. The droplet diameter is 0.4 and it oscillates in a 2 × 2
square domain. The initial shape of the droplet is given in polar coordinates (r, θ ) by
r (θ) = r0 + cos(nθ),
where n is the integer mode of oscillation. The theoretical small-amplitude inviscid oscillation
frequency is then the following (Lamb 1971):
ωn2 =
(n 3 − n)σ
(ρd + ρe )r03
where ρd and ρe are the droplet and exterior densities. This test case is similar to that
used by Torres and Brackbill (2000) and Hermann (2008) and we compare the results
obtained using Gerris with the results presented in both of these papers. A comparison
between the theoretical and numerical values of the frequency is given in figure 2 for n = 2
and = 0.05.
Note that the results of Torres and Brackbill and Herrmann were obtained using a small
explicit viscosity and for density ratios of 1/100. This test case was initially selected by Torres
and Brackbill to demonstrate the importance of accurate surface tension representation. Given
the low explicit viscosity, any spurious currents due to inaccurate surface tension will not
be damped and will strongly affect the solution. The accuracy of Gerris in this respect is
confirmed by this test case. Another measure of the accuracy of the solution is the estimation
of the numerical damping induced by the various discretization errors of the numerical
Fluid Dyn. Res. 41 (2009) 065001
Equivalent Laplace number
Kinetic energy
D Fuster et al
100 000
10 000
Number of grid points per droplet diameter
Figure 3. (a) Evolution of the total kinetic energy for an oscillating droplet for the spatial
resolutions given in the legend (in grid points per diametre). (b) Equivalent Laplace numbers.
scheme. The amount of numerical damping can be estimated by computing an equivalent
viscosity. With viscosity, kinetic energy is expected to decrease as
where C is a dimensionless constant (set empirically to 30), ν the viscosity, D the
droplet diameter and t the time. Using curve fitting the damping coefficient b = Cν/D 2
can be estimated (solid lines in figure 3(a)). An equivalent Laplace number can then be
computed as
La ≡
σ C2
= 2 3.
ρb D
The numerical viscosity and thus the equivalent Laplace number depend on spatial resolution
as illustrated in figure 3(b). The kinetic energy is seen to behave in a physically consistent
manner at all resolutions. High Laplace numbers can be reached at high resolution, which
confirms that the overall stability of the method is not due to an excess of numerical viscosity.
In what follows we demonstrate how the combination of adaptivity, accurate surface
tension and overall robustness of the method leads to novel applications by focusing on
three problems: 2D atomization, retraction of liquid sheets and small-scale air–water wave
3. Atomization
Atomization processes are of significant interest in various industrial and research fields.
However, our understanding of the physical processes leading to the break-up of liquid masses
is still far from complete. Although different phenomena have been separately studied in the
literature, the complexity of the analysis of the whole processes has prevented investigators
from obtaining a complete understanding of all the processes underlying atomization.
Nowadays, numerical simulations have become a major avenue of investigation
for atomization processes, alongside theory and experimentation. Thanks to the recent
developments introduced earlier in this paper, it is now possible to perform simulations with
a high level of accuracy in conditions very close to those of atomization experiments.
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Table 2. Simulation conditions for the analysis of atomization processes.
U (m s−1 )
ρ (kg m−3 )
µ (Pa s)
σ (N m−1 )
1.7 × 10−5
Generally, three different parts of the processes have been investigated separately, namely,
the flow inside the injector, the primary atomization zone near the injector and the secondary
atomization zone further downstream. The flow inside the injector establishes the input
boundary conditions at the chamber. In the primary atomization zone, the liquid jet that is
injected into the chamber is exposed to the high-velocity gas flow. The perturbations induced
at the interface quickly grow downstream, becoming nonlinear and producing long ligaments
that eventually break into droplets. Finally, in the secondary atomization zone, large fragments
of the central jet zone and droplets created upstream in the primary atomization region may
be further broken into droplets or coalesce and the undisturbed jet core can no longer be
Different spatial scales appear across the domain. In some parts, such as near the injector,
a fine mesh resolution is required to capture the correct wave amplification, the ligament
formation or the ligament break-up. In other parts, such as in the secondary atomization zone,
fine mesh resolution is required to capture the whole range of droplet sizes generated there
or carried from upstream. However, with current computational resources, such a degree of
resolution cannot be applied to the entire domain and AMR techniques are required in order
to concentrate the computational effort just in those zones where it is needed.
When the whole process is investigated, one of the main issues is the presence of strong
interactions between different parts of the system: the flow upstream, which is inside the
injector nozzle, determines the parameters of the atomization at the entrance of the chamber,
such as the intensity of turbulence introduced in the system or the size of the boundary layers.
The flow in the secondary atomization zone can create some vortices and droplets. These can
be sent back upstream by recirculating flow and influence the flow in the primary atomization
zone. The present work aims at capturing all the interactions between the different parts of
the system. Although some assumptions and simplifications are still required, it is shown
that current advances in the numerical simulation of atomization processes are enough to
capture important information and to recover quantitative experimental measurements for real
air/water systems.
3.1. Simulation of atomization processes in conditions close to experiments
In this section, the effect of the injector in the primary atomization zone is investigated.
In order to perform simulations with large resolution in big domains, the 2D problem is
3.1.1. Simulation setup. Conditions contained in table 2 are used in this work to show the
potential of the code to handle air/water simulations.
These conditions, which define the characteristic dimensionless parameters contained in
table 3, correspond to the experimental investigation of air/water atomization processes at
large momentum ratios M = ρl Ul2 /(ρg Ug2 ) reported by Ben Rayana (2007).
The simulation setup is depicted in figure 4. Although the main goal is to analyze the
characteristic frequencies and flow patterns near the injector, it has been found that both
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Table 3. Representative dimensionless numbers for the conditions given in table 2, where µ is
the viscosity, δ is the thickness of the boundary layer at the entrance, ρ is the density and σ is the
surface tension.
W e∗g
µg /µl
ρg /ρl
ρl Ul Hl /µl
ρg Ug δg /µg
ρg (Ug )2 δg /σ
1.2 × 10−3
ρg Ug2
ρl Ul2
Figure 4. Planar liquid sheet (bottom) atomized by an upper high-velocity gas jet when a
statistically steady state is reached. The thickness of each of the jets is 1 cm. Colors represent
the vorticity field intensity in the counterclockwise direction (red) and the clockwise direction
(blue). Domain size: 11 cm × 6 cm. The origin of the x-coordinate is placed at the border of the
injector, where the two layers are in contact. The origin of the y-coordinate is placed at the bottom
of the domain.
parameters can be much influenced by the domain boundaries if the domain is not large
enough. If the output boundary condition is imposed before the jet is definitely broken into
droplets, some wave reflection can be induced at the outflow boundaries and they significantly
modify the characteristic frequencies observed some distance downstream. In order to avoid
this effect, a large domain of 16 × 6 cm is considered. The liquid and gas heights of the sheet
at the entrance are, for both phases, Hg = Hl = 1 cm. At the entrance, two boundary layers are
introduced to reproduce the velocity profiles near the solid walls (δl = δg = 179 µm). Between
the liquid and gas jets, a separator plate with a minimum thickness of 150 µm and a slope (or
the opening half-angle) of β = 7◦ is included where non-slip boundary conditions at the walls
are imposed (the separator plate can be seen more clearly in the zoomed view of figure 7).
The inclusion of the separator plate inside the simulation domain is important to capture some
of the effects of the injector geometry on the flow (e.g. to capture the effect of the thickness
of the separator plate e and the opening half-angle β).
Finally, a standard outflow boundary condition (zero normal derivative of normal velocity
and zero pressure) is imposed at the right boundary of the domain, whereas slip boundary
conditions are set at the top and bottom boundaries as well as on the left boundary in those
zones where there is no entering jet.
Regarding the mesh, an automatic octree refinement based on two different criteria has
been used: the gradient of the VOF color function and the norm of the vorticity. The first
criterion allows us to control the degree of accuracy with which the interface is captured. The
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Figure 5. Simulation domain and one example of a typical mesh distribution.
second one controls the level of resolution achieved to capture the turbulent vortices. Both
criteria are also space-dependent in order to concentrate the computational effort on the zone
near the injector. In particular, the maximum spatial resolution in the different parts of the
domain is
9.76 µm for x < 2 cm and y < 3 cm,
39 µm for 2 < x < 3 cm and y < 3 cm,
156 µm for x > 3 cm and y < 3 cm.
For y > 3 cm, the interface is refined with cells of 156 µm if some liquid reaches this
zone. If no liquid is present, a coarse mesh is used to dampen the vortices entering
this zone.
In addition, the grid near the solid walls is maximally refined (9.76 µm) in order to maximize
the accuracy of the simulation near the separator and to capture the generation of instabilities
by the trailing edge of the separator plate.
We would like to stress the importance of octree refinement (figure 5).
• On the one hand, it has allowed us to consider a domain where the boundaries are
relatively far from the injector without a prohibitive increase of the number of cells. For
distances further than 2 cm from the injector, the mesh is progressively coarsened. Thus,
although the accuracy of the results in this zone is not as good as near the injector, the
influence of the boundary conditions on the observed frequencies is negligible.
• On the other hand, very refined meshes can be used near the injector (of the order of
10 µm), which allows us to perform high-accuracy simulations of the flow in this zone.
3.1.2. Analysis of the simulation results. The new simulations allow us to capture various
flow patterns in the different parts of the domain (figure 4). At the trailing edge of the separator
plate, the liquid is exposed to the high-velocity gas sheet and the exchange of momentum
between the gas and liquid jets accelerates the liquid. The amplitude of the waves generated
at the interface is significant even immediately downstream of the trailing edge.
Qualitatively, we can observe the amplification of the instabilities induced near the
injector as well as the generation of ligaments, which are finally broken into droplets. The
presence of the liquid potential cone is also clearly distinguished. In figure 6, the superposition
of the graph of the interface during 0.2 s reveals that the length of the potential cone is between
3 and 4 cm, which agrees with the experimental correlation found by Ben Rayana (2007):
L cone = 12M −1/2 HL = 3 cm.
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Figure 6. Dots: projection of the position of the interface during 0.2 s. Line: slope of the potential
cone extracted from the experimental correlation.
Figure 7. Zoom of the flow near the separator plate. Sequence of images of the generation, growth
and propagation of the perturbation induced near the injector. Time lag between images: 5 µs.
yr represents the distance from the separator plate in the y-direction in cm.
Focusing the analysis on the zone near the injector, we can observe waves of relatively
large amplitude immediately behind the injector (of the order of the boundary layer thickness
and the separator plate thickness). The generation of these waves turns out to be a periodic
process summarized in the sequence of images depicted in figure 7 and described below:
After one wave has been advected downstream, the liquid level in the injector
surroundings is approximately constant and it is sustained for a relatively long distance
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Distance from the separator plate: 0.4 mm
yr (mm)
yr (mm)
Distance from the separator plate: 0.2 mm
U (m s–1)
U (m s–1)
Figure 8. Two examples of the instantaneous profiles of the x velocity component at two different
distances from the trailing edge of the separator plate when the interface is at the level of
the injector (+ symbols) and when the liquid level has decreased due to the amplification and
propagation of a large wave downstream (× symbols).
compared with the thickness of the separator plate. At this moment, some perturbations are
induced at the interface due to the turbulence generated in the gas phase near the interface.
These disturbances can have an amplitude as large as the thickness of the separator plate.
When a disturbance occurs, the wave amplifies and once the crest of the wave goes beyond
the gas boundary layer, a sudden acceleration and amplification is produced. Due to the large
amplitude of the wave propagating downstream, the mass conservation makes the liquid level
around the injector decrease. The volume that was initially occupied by the liquid is now
filled with a recirculating gas flow. During these moments, the effective thickness of the gas
boundary layer is much larger than the amplitude of the disturbances induced at the interface
(figure 8). The decrease of the liquid level is large enough to avoid any wave reaching the gas
bulk. In addition, the large gas vortices produced in the gas boundary layer make the waves
propagate in both directions, upstream and downstream. At this moment, the disturbances
induced at the interface are not in contact with the gas jet and are not amplified. It is not until
the liquid level near the injector is recovered that the next wave of large amplitude appears.
When a big wave is generated, some ligaments appear that are eventually broken into
relatively large droplets (see figure 9). This process governs the appearance of droplets near
the injector and will be discussed in more detail in the next section.
Then, we can conclude that the current simulations have been shown to capture a realistic
behavior of the atomization processes for air/water systems at large momentum ratios. This
result is remarkable if we take into account the shortage of numerical simulations for this kind
of system due to the different numerical issues related with the large density and viscosity
ratios encountered at the interface. The presence of the separator plate has made it possible to
introduce correct boundary conditions at the entrance, and thanks to AMR, large domains can
be simulated, making it possible to diminish the influence of the boundary conditions on the
dynamics near the separator plate.
4. Capillary retraction of a liquid sheet
4.1. Introduction
During the atomization, drops may be formed by several distinct mechanisms. A general
understanding of these processes is still lacking and is at the heart of many fundamental
studies on atomization: the Kelvin–Helmholtz instabilities, jet and filament break-up, and
liquid sheet deformation. In particular, the destabilization of a liquid sheet is known to detach
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Figure 9. Generation of a single drop from a ligament in the primary atomization zone. For
the sake of clarity, the velocity vectors are just depicted at level of refinement 8 (78 µm).
Top: global view. Bottom: a zoom near the topmost part of the interface.
small droplets in situations such as the crown splash illustrated by the iconic photographs of
Edgerton (1987).
A particular feature of liquid sheet dynamics is that its free ends retract because of
the surface tension, thus forming a rim that usually grows with time and is often unstable.
As explained and shown by Taylor (1959) and Culick (1960), the retraction dynamics
of a homogeneous liquid sheet initially at rest reaches a constant speed, the so-called
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Table 4. Characteristic parameters used to write the rest of the variables in dimensionless form:
characteristic length L c , characteristic velocity Uc and characteristic density ρc .
L c (m)
Uc (m s−1 )
ρc (kg m−3 )
5 × 10−3
Taylor–Culick velocity vTC = 2σ/(ρe), where σ is the surface tension, ρ the liquid density
and e the thickness of the liquid sheet. This velocity comes from the momentum balance
between the liquid entering the rim and the surface tension. A full description of the dynamics
of the liquid sheet requires 3D analysis and extensive numerical simulations and is beyond
the scope of the present paper. Here, we will instead focus on the 2D retraction dynamics.
Although we know that the 2D flows thus obtained may be unstable in 3D (Fullana and
Zaleski 1999), the properties of the 2D retraction are remarkable and a better understanding
of these dynamics is needed to improve our understanding of rim instability even in 3D.
Rich and complex behaviors have already been observed in 2D, such as capillary waves
and a pinching of the film that becomes more pronounced with surface tension (Song and
Tryggvason 1999) or viscous-capillary-balance retraction speed (Debrégeas et al 1995, 1998,
Brenner and Gueyffier 1999). Here, we will show in particular that the 2D film retraction does
not pinch the liquid film and that it becomes unstable at low Ohnesorge numbers.
4.2. Numerical simulation
Our aim is thus to numerically investigate the motion and stability of the edge of a liquid sheet
in 2D. We consider a 2D flow and we try to be as close as possible to the typical flows of an
air/water system. The capillary retraction is controlled by the density and viscosity ratios,
and by the Ohnesorge number (Oh =√µ/(ρσ l0 )1/2 ), where l0 is the initial homogeneous
thickness of the sheet. Note that Oh = lvc /l0 where lvc = µ2 /(ρσ ) is the viscous-capillary
length. For water, note that lvc ' 14 × 10−9 m, a particularly small length. For a sheet with
l0 ' 50 × 10−6 m (five times the finest grid size used in the simulation of section 3.1.2),
we thus have Oh ' 0.017. We perform the study in the reference frame where the film edge
is at rest to allow large computation times. In this framework, both the liquid and air are
entering into the calculation domain (at the left boundary, see figure 10) at the Taylor–Culick
velocity vTC .
The dynamical equations are solved using the dimensionless version of the equations.
The physical parameters used in the simulations to obtain the dimensionless system are
summarized in table 4.
We use a computational domain of size 4 × 1 in dimensionless units. The initial mesh has
4 × 27 × 27 cells. The mesh is adapted based on two criteria, namely the interface position and
the magnitude of the vorticity field, in order to maintain a good resolution around the interface
as well as wherever the vorticity is large.
Finally, a Dirichlet boundary condition is imposed on the velocity field on the left
boundary where the gas and liquid are penetrating while an outlet condition is taken on the
right-hand side. The density ratio is 850 and the viscosity ratio is 50, which corresponds to
an air/water system. In our dimensionless units, the initial length of the sheet is 18e0 and
the unit-square subdomain corresponds to 20e0 , where e0 = l0 /L c is the initial dimensionless
thickness of the sheet (see figure 10). We introduce the liquid capillary time τ =
which allows us to define the dimensionless time t ∗ = τt .
ρ L e03 /σ ,
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Figure 10. Spatial domain of the retracting sheet. The domain consists of 4 × 1 square
subdomains. The liquid sheet is surrounded by a gas with a density 850 times smaller.
Figure 11. Evolution of the sheet interface at t ∗ = 0.00, 5.29, 10.58, 15.87, 21.17, 31.75, 42.33
and 52.92 (from top to bottom) for four different Ohnesorge numbers: 0.1 (a), 0.055 (b), 0.01 (c)
and 0.005 (d). A rim is formed at the edge of the free end of the sheet and a neck appears just
behind the rim.
4.3. Results
Figure 11 shows snapshots of the liquid sheet evolution at eight different times and for four
different Ohnesorge numbers ranging from 0.005 to 0.1. We observe that after a short transient
(a few units of times), the rims almost stay centered on the same position as time evolves.
This demonstrates that the rim recedes at the Taylor–Culick speed vTC . A neck, more or less
pronounced, forms behind the rim.
This figure illustrates also an instability mechanism in the retraction dynamics for small
Ohnesorge numbers. Indeed, for high Ohnesorge numbers (0.1 and 0.055 in figure 11), viscous
effects are important and the rim is stable (see figures 11(a) and (b)). On the other hand, as
the Ohnesorge number decreases, the evolution of the rim is significantly different. As shown
in figures 11(c) and (d), the rim is no longer circular but elongated in the x-direction and
oscillations of the rim are observed.
Moreover, the neck is followed in the sheet by a stationary capillary wave whose
amplitude is exponentially decreasing in space for Oh 6 0.5, which is always the case in
our simulations (Agbaglah et al, in preparation; Gordillo, private communication 2007). The
Fluid Dyn. Res. 41 (2009) 065001
Oh = 0.1
Oh = 0.01
Oh = 0.005
Oh = 0.1
Oh = 0.01
Oh = 0.005
Rim thickness/e0
D Fuster et al
Figure 12. (a) Evolution of the size of the rim for Oh = 0.1, Oh = 0.01 and Oh = 0.005. (b)
Evolution of the neck thickness for Oh = 0.1, Oh = 0.01 and Oh = 0.005. The sheet is stable for
Oh = 0.1 and becomes unstable for Oh smaller than 0.01.
stationary capillary waves structure along the liquid sheet is more pronounced for low
Ohnesorge numbers as can be observed in figure 11, but it is important to note that no breakup of the liquid sheet has been observed in the numerics, in contrast to the prediction of
finite-time singular pinch-off (Song and Tryggvason 1999). The appearance of a finite-time
singularity is well known for 3D capillary dynamics such as liquid jet pinch-off for instance
(Eggers 1997). Here, however, the surface tension cannot enhance the finite time singularity in
2D and the neck is dominated by the balance between capillarity, inertia and viscous stress.
Figure 12, where the rim and neck thicknesses are plotted for several Ohnesorge numbers,
also illustrates the change of behavior in the retraction dynamics. While the large Ohnesorge
number (0.1 here) shows the expected square-root-of-time regime (Keller et al 1995), the
elongation of the rim is obvious for smaller Oh numbers (0.01 and 0.005 in figure 11).
Similarly, figure 12(b) shows a smooth convergence of the neck thickness for the highest
Ohnesorge number while strong and irregular oscillations of this thickness can be observed
for the smaller Ohnesorge numbers. This instability is also illustrated in figure 13, where
the neck thikness is shown as a function of time for Ohnesorge numbers close to the
instability threshold.
Finally, figure 14 shows that the instability of the edge of the sheet is in fact related to the
presence of vortices nucleated by the dynamics. For low Ohnesorge numbers, the vorticity,
created by the neck and capillary wave structures, is advected inside the rim and makes it
oscillate. Such highly vortical flows can indeed destabilize the rim (see Song and Tryggvason
1999). For large Ohnesorge numbers, the vorticity remains concentrated around the neck and
no instability can develop. This suggests that the instability can be seen as a particular case of a
Kelvin–Helmholtz instability, due to the shear between the gas and the retracting liquid sheet.
Indeed, the oscillation frequency observed in both figures 12 and 13 corresponding to ω ∼ 1
is in good agreement with the Kelvin–Helmholtz mode calculated with the Taylor–Culick
retraction velocity and the wavelength of spatial oscillation of the film.
From our set of numerical simulations, we also deduce that this instability develops for
Ohnesorge numbers smaller than 0.013 ± 0.003 (figure 13). For an air–water system, this
would correspond to a film thickness approximately larger than 100 µm.
5. Wave breaking
Aside from the combustion-engine injectors discussed earlier, atomization of liquid–gas
interfaces is also typical of the breaking of steep air–water waves. This is an important
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Oh = 0.02125
Oh = 0.015625
Oh = 0.01
Figure 13. Evolution of the neck thickness for different Oh values. A transition regime is
encountered for Oh between 0.016 and 0.01.
problem in geophysics since the dynamics of water waves controls a number of fundamental
processes including the transfer of soluble gases between the atmosphere and the ocean and
the formation of atmospheric aerosols. Depending on wave steepness, various regimes can be
distinguished: linear waves, weakly nonlinear waves, stable nonlinear waves, spilling waves,
plunging breakers, etc. Also—while the water phase largely controls the transfer of slightly
soluble gases for weakly nonlinear waves—the air phase must be taken into account when
dealing with steep spilling or breaking waves. Small-scale waves are predominant at the ocean
surface and may be responsible for an important fraction of the global ocean–atmosphere
fluxes. At these scales, surface tension plays an important role.
From a numerical analysis perspective, two main trends can be distinguished:
1. Potential flow theory in conjunction with boundary-integral methods has been
successfully applied to the study of the long-time evolution of large-scale water waves
and their interactions (Dias and Bridges 2006).
2. Direct numerical simulation (DNS) techniques initially developed in the context of
channel-flow turbulence studies have been extended to flows beneath a fixed or weakly
deformed free surface (Tsai et al 2005).
Both of these approaches have important limitations. Potential flow techniques cannot
describe the boundary layers controlling gas exchange while existing DNS methods do not
include the nonlinear deformations of the interface and/or are limited to a single phase
(Tsai et al 2005, Tsai and Hung 2007). Further advances could be made by using recently
developed methods for the direct simulation of two-phase interfacial flows. While a few
studies using these techniques have been published, they are often restricted in scope
due to two main limitations: inaccuracies of the numerical methods (particularly regarding
momentum conservation and surface tension representation) and prohibitive computational
cost at high spatial resolutions (particularly in 3D). Both these limitations are addressed in
the recent developments discussed earlier. This section illustrates the potential of these new
techniques for the study of steep nonlinear air–water waves.
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Figure 14. Vorticity field at t ∗ = 5.29 for different Ohnesorge numbers (from top to bottom: 0.1,
0.055, 0.01 and 0.005). A pair of vortices appears inside the rim only for Oh = 0.01 and 0.005.
An example of the results obtained for the evolution of
√ a small-amplitude viscous Stokes
wave is given in figure 15. Non-dimensional times t 0 = t g/λ with λ being the wavelength
are 1.32, 2.65, 3.97 and 5.30 from top to bottom. Surface tension is included, which leads to
the formation of capillary waves riding in front of the gravity wave. No external forcing is
added and the amplitude of the wave slowly decreases due to viscous dissipation. Both the
thin viscous boundary layer in the water phase and the thicker boundary layer in the gas phase
are evident with strong vorticity generation associated with the high curvatures of capillary
waves. These results are similar (for the water phase) to those obtained by Tsai and Hung
(2007) using their single-phase free-surface solver.
As with other implicit front-capturing methods, the interface topology is not restricted:
air entrainment and spray formation can thus be modeled directly provided resolution is high
enough. Figure 16 illustrates preliminary results for the simulation of a 3D unstable steep
Stokes wave. To illustrate 3D breaking, the initial third-order Stokes wave steepness parameter
ak is varied from 0.35 to 0.45 in the direction transverse to the wave motion. The boundary
conditions are periodic in the direction of wave motion and symmetrical in the transverse
direction. A characteristic plunging breaker is formed at t 0 = 2, and then leads to the trapping
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Figure 15. p
Evolution of a small-amplitude third-order Stokes wave with surface tension. ak =
0.28, Re = gλ3 /ν = 40 000, Bo = ρλ2 g/σ = 333, air/water density and viscosity ratios. ρ and
ν are the water density and viscosity and σ is the water/air√surface tension coefficient. Color
represents vorticity (dark blue is −50; dark red is +50 times g/λ). The wave is traveling from
left to right.
of a large air ‘tube’ and the formation of a secondary plunging sheet at t 0 = 2.5. The edge of
this sheet is rounded due to capillary retraction. Eventually, holes form immediately behind
the sheet’s edge and grow due to capillary retraction, leading to the formation of cylindrical
‘fingers’, which are themselves unstable due to Plateau–Rayleigh instability. At t 0 = 2.8 the
sheet has impacted with the front of the wave and further atomization occurs as the droplets
themselves impact on the interface.
As in the previous simulation, the octree mesh is adapted at every time step according to
the local magnitude of vorticity and interface curvature. Figure 17 illustrates a detail of the
variable-resolution VOF representation of the interface. The maximum resolution allowed is
set to 1/λ = 1/1024. For the Reynolds and Bond numbers of 40 000 and 33 333 chosen,
this guarantees that both the smallest viscous and capillary scales are resolved. This is
confirmed by the observation (figure 17 for example) that most droplets and bubbles have
sizes significantly larger than the mesh size.
Adaptivity is crucial for obtaining the results presented here. Figure 18 illustrates the
evolution with time of the total number of grid points. The initial wave is smooth and can
be resolved accurately with a limited number of elements. As the plunger forms, the local
interface curvature increases and requires higher resolution. Plunger impact at t 0 = 2.1 creates
other highly curved structures and, finally, sheet atomization at t 0 = 2.5 further increases the
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
t = 2
t = 2.5
t = 2.8
Detail at t = 2.8
Figure 16. Evolution of the plunging breaker created by the instability of a steep Stokes
wave. 0.35 < ak < 0.45, Re = 40 000, Bo = 33 333, air/water density
√ and viscosity ratios. Color
represents the magnitude of velocity (dark blue is zero; dark red is gλ).
resolution requirements. The final number of grid points of ≈1.4 million can be compared to
the ≈1 billion grid points that would be required if a regular Cartesian grid was used. Even a
method using adaptivity but with a constant interface resolution would require about 50 times
more grid points.
6. Conclusion
Current advances in numerical simulation are allowing us to go deep into the understanding of
the physical processes underlying multiphase flows. The present work presents three related
examples showing the potential of current methods: the investigation of the instabilities in
atomization processes, the capillary retraction of a liquid sheet and 2D and 3D wave breaking.
The accurate numerical schemes described above together with AMR techniques have
made such simulations possible even for the difficult regimes of air/water systems.
Regarding atomization, a liquid sheet subjected to the effect of a high-velocity jet is
considered. The concentration of the mesh cells in the zone near the injector has allowed us
to accurately capture the flow structures in this zone, whereas further away, the mesh has
been coarsened in order to reduce the computational effort. Together with the removal of the
small droplets introduced in the system, this has allowed us to perform accurate simulations in
conditions used experimentally to investigate air/water atomization processes. Quantitatively
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Figure 17. Detail of the adaptive VOF interface representation for the breaking wave at t 0 = 2.8.
Number of grid points
900 000
700 000
Plunger impact
500 000
300 000
Plunger formation
100 000
Figure 18. Evolution with time of the total number of grid points for the 3D wave breaking
the length of the potential cone correctly fits the experimental measurements. In addition,
new mechanisms leading to the appearance of characteristic frequencies have been captured.
The new simulations seem to point out that the appearance of the large waves induced in the
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
system are related to the characteristic time required by the flow near the injector to recover
its level once the previous wave is propagated downstream. This regime is probably linked to
a high momentum ratio M or to high ratios Ug /Ul .
As a second example, the 2D retraction of a water sheet in air is considered. Again,
even for large density and viscosity ratios, no numerical instabilities are found. The new
simulations have shed some light into the drop-formation process. Two distinct regimes are
found for the neck contraction process as a function of the Ohnesorge number. For values
larger than 0.016, a smooth neck contraction is observed and its evolution tends to a constant
value. For values smaller than 0.01, a turbulent flow is predicted inside the rim, which
produces large oscillations in both the neck thickness and the rim shape.
Finally, the preliminary study of air/water wave breaking further illustrates the
robustness, accuracy and efficiency of the numerical methods presented. In this case,
adaptivity is shown to provide a gain of three orders of magnitude in simulation size compared
to an equivalent regular-Cartesian-grid simulation. Allowing adaptivity along the interface
based on interface curvature reduced the simulation size by a factor of 50 compared to an
adaptive technique resolving the interface at constant resolution.
Almgren A S, Bell J B, Colella P, Howell L H and Welcome M L 1998 A conservative adaptive projection method
for the variable density incompressible Navier–Stokes equations J. Comput. Phys. 142 1–46
Aulisa E, Manservisi S, Scardovelli R and Zaleski S 2007 Interface reconstruction with least-squares fit and split
advection in three-dimensional cartesian geometry J. Comput. Phys. 225 2301–19
Ben Rayana F 2007 Étude des instabilités interfaciales liquide-gaz en atomisation assistée et tailles de gouttes PhD
Thesis Institut National Polytechnique de Grenoble
Brenner M and Gueyffier D 1999 On the bursting of viscous films Phys. Fluids 11 737–9
Chorin A J 1969 On the convergence of discrete approximations to the Navier–Stokes equations Math. Comput. 23
Coantic M 1980 Mass transfert across the ocean–air interface: small scale hydrodynamic and aerodynamic
mechanisms PhysicoChem. Hydrodyn. 1 249–79
Compère G, Marchandise E and Remacle J-F 2008 Transient adaptivity applied to two-phase incompressible flows
J. Comput. Phys. 227 1923–42
Culick F E C 1960 Comments on a ruptured soap film J. Appl. Phys. 31 1128–9
Dai M and Schmidt D P 2005 Adaptive tetrahedral meshing in free-surface flow J. Comput. Phys. 208 228–52
Debrégeas G, de Gennes P and Brochard-Wyart F 1998 The life and death of bare viscous bubbles Science 279
Debrégeas G, Martin P and Brochard-Wyart F 1995 Viscous bursting of suspended films Phys. Rev. Lett. 75 3886–9
Dias F and Bridges T J 2006 The numerical computation of freely propagating time-dependent irrotational water
waves Fluid Dyn. Res. 38 803–30
Duchemin L, Popinet S, Josserand C and Zaleski S 2002 Jet formation in bubbles bursting at a free surface Phys.
Fluids 14 3000
Edgerton H 1987 Stopping Time (New York: Harry N Abrams)
Eggers J 1997 Nonlinear dynamics and breakup of free-surface flows Rev. Mod. Phys. 69 865–930
Eggers J and Villermaux E 2008 Physics of liquid jets Rep. Prog. Phys. 71 036601
François M M, Cummins S J, Dendy E D, Kothe D B, Sicilian J M and Williams M W 2006 A balanced-force
algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework
J. Comput. Phys. 213 141–73
Fullana J M and Zaleski S 1999 Stability of a growing end-rim in a liquid sheet of uniform thickness Phys. Fluids 11
Fuster D, Bagué A, Boeck T, Le Moyne L, Leboissetier A, Popinet S, Ray P, Scardovelli R and Zaleski S 2009
Simulation of primary atomization with an octree adaptive mesh refinement and VOF method Int. J. Multiphase
Flow 35 550–65
Fyfe D E, Oran E S and Fritts M J 1988 Surface tension and viscosity with Lagrangian hydrodynamics on a triangular
mesh J. Comput. Phys. 76 349–84
Greaves D M 2004 A quadtree adaptive method for simulating fluid flows with moving interfaces J. Comput. Phys.
194 35–56
Fluid Dyn. Res. 41 (2009) 065001
D Fuster et al
Gueyffier D, Nadim A, Li J, Scardovelli R and Zaleski S 1998 Volume of fluid interface tracking with smoothed
surface stress methods for three-dimensional flows J. Comput. Phys. 152 423–56
Herrmann M 2008 A balanced force refined level set grid method for two-phase flows on unstructured flow solver
grids J. Comput. Phys. 227 2674–706
Josserand C and Zaleski S 2003 Droplet splashing on a thin liquid film Phys. Fluids 15 1650–7
Kang M, Fedkiw R P and Liu X D 2000 A boundary condition capturing method for multiphase incompressible flow
J. Sci. Comput. 15 323–60
Keller J, King A and Ting L 1995 Blob formation Phys. Fluids 7 226–8
Lamb H 1971 Hydrodynamics (Cambridge: Cambridge University Press)
Popinet S 2003 Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries
J. Comput. Phys. 190 572–600
Popinet S 2007 The Gerris Flow Solver http://gfs.sf.net
Popinet S 2009 An accurate adaptive solver for surface-tension-driven interfacial flows J. Comput. Phys. 228
Popinet S and Zaleski S 1999 A front-tracking algorithm for accurate representation of surface tension Int. J. Numer.
Methods Fluids 30 775–93
Renardy Y and Renardy M 2002 PROST—A parabolic reconstruction of surface tension for the volume-of-fluid
method J. Comput. Phys. 183 400–21
Rudman M 1998 A volume-tracking method for incompressible multifluid flows with large density variations Int. J.
Numer. Methods Fluids 28 357–78
Scardovelli R and Zaleski S 1999 Direct numerical simulation of free-surface and interfacial flow Annu. Rev. Fluid
Mech. 31 567–03
Shin S, Abdel-Khalik S I, Daru V and Juric D 2005 Accurate representation of surface tension using the level contour
reconstruction method J. Comput. Phys. 203 493–516
Song M and Tryggvason G 1999 The formation of thick borders on an initially stationary fluid sheet Phys. Fluids 11
Sussman M, Smereka P and Osher S 1994 A level set approach for computing solutions to incompressible two-phase
flow J. Comput. Phys. 114 146–59
Sussman M, Smith K M, Hussaini M Y, Ohta M and Zhi-Wei R 2007 A sharp interface method for incompressible
two-phase flows J. Comput. Phys. 221 469–505
Taylor G I 1959 The dynamics of thin sheets of fluid III. Disintegration of fluid sheets Proc. R. Soc. A 253 313–21
Torres D J and Brackbill J U 2000 The point-set method: front-tracking without connectivity J. Comput. Phys. 165
Tsai W T, Chen S M and Moeng C H 2005 A numerical study on the evolution and structure of a stress-driven
free-surface turbulent shear flow J. Fluid Mech. 545 163–92
Tsai W and Hung L 2007 Three-dimensional modeling of small-scale processes in the upper boundary layer bounded
by a dynamic ocean surface J. Geophys. Res. 112 C02019
Yarin A L and Weiss D A 1995 Impact of drops on solid surfaces: self-similar capillary waves and splashing as a new
type of kinematic discontinuity J. Fluid Mech. 283 141–73
Fly UP