# A momentum-conserving, consistent, Volume-of-Fluid method for incompressible
flow on staggered grids

D. Fuster  T. Arrufat  M. Crialesi-Esposito  Y. Ling  L. Malan  S. Pal  R.
Scardovelli  G. Tryggvason  S. Zaleski  Sorbonne Université et CNRS,
Institut Jean Le Rond d’Alembert, UMR 7190, Paris, France  CMT-Motores
Térmicos, Universitat Politécnica de Valéncia, Camino de Vera, s/n, Edificio
6D, Valencia, Spain  InCFD, Dept. of Mechanical Engineering, University of
Cape Town, South Africa  DIN - Lab. di Montecuccolino, Università di Bologna,
I-40136 Bologna, Italy  Mechanical Engineering, Johns Hopkins University,
Baltimore, USA  Dept. of Mechanical Engineering, Baylor University, Waco, TX,
USA

###### Abstract

The computation of flows with large density contrasts is notoriously
difficult. To alleviate the difficulty we consider a consistent mass and
momentum-conserving discretization of the Navier-Stokes equation.
Incompressible flow with capillary forces is modelled and the discretization
is performed on a staggered grid of Marker and Cell type. The Volume-of-Fluid
method is used to track the interface and a Height-Function method is used to
compute surface tension. The advection of the volume fraction is performed
using either the Lagrangian-Explicit / CIAM (Calcul d’Interface Affine par
Morceaux) method or the Weymouth and Yue (WY) Eulerian-Implicit method. The WY
method conserves fluid mass to machine accuracy provided incompressiblity is
satisfied which leads to a method that is both momentum and mass-conserving.
To improve the stability of these methods momentum fluxes are advected in a
manner “consistent” with the volume-fraction fluxes, that is a discontinuity
of the momentum is advected at the same speed as a discontinuity of the
density. To find the density on the staggered cells on which the velocity is
centered, an auxiliary reconstruction of the density is performed. The method
is tested for a droplet without surface tension in uniform flow, for a droplet
suddenly accelerated in a carrying gas at rest at very large density ratio
without viscosity or surface tension, for the Kelvin-Helmholtz instability,
for a falling raindrop and for an atomizing flow in air-water conditions.

###### keywords:

Multiphase Flows , Navier-Stokes Equations , Volume of Fluid , Surface Tension
, Large Density Contrast

††journal: Computers and Fluids

##  1 Introduction

Multiphase flows abound in nature, but their stable and accurate computation
remains elusive in many cases. As a case in point, many numerical methods used
for two-phase incompressible flow are strongly unstable for large density
contrasts and large Reynolds numbers. Experience with such simulations shows
that the presence of surface tension is an aggravating factor. The large
density contrasts that are of interest are air/water or gas/liquid-metal, with
$\rho_{l}/\rho_{g}$ of the order of $10^{3}$ or $10^{4}$. The large density
contrasts are a difficulty whether one deals with any of the three major
interface advection methods, Level-Set, Volume-of-Fluid (VOF) or Front-
Tracking, or with combinations such as CLSVOF. (The term density contrast is
perferable to density ratio since it encompasses ratios both much larger than
one and much smaller than one.)

Several methods have been used to alleviate the high-density-contrast
difficulties. It has been observed by several authors that making the
momentum-advection method conservative improves the situation. For
incompressible flow, momentum-conserving methods have been initially proposed
by [1], and by several other authors since [2, 3, 4, 5, 6]. These methods have
been shown to improve the stability of the numerical results in various
situations. In particular, liquid-gas flows with very contrasted densities, as
for example in the the process of atomization, cause serious problems that are
resolved by using momentum-conserving methods. In that case another oft-
suggested solution is to increase the number of equations from the standard
four equations to five, six or seven equations, by introducing new field
variables in each phase. The addition of one more density $\rho_{i}$, momentum
$\rho_{i}{\bf u}_{i}$ or energy variable $\rho_{i}e_{i}$ increases the number
of equations. The authors of ref. [7] used seven equations, those of [8] used
five equations and six equations were used in [9]. The last three references
also use a momentum-conserving formulation. Several authors, including some of
those cited above, have argued that the difficulty may come from gas
velocities of order $u_{g}$ being mixed with liquid densities of order
$\rho_{l}$. Both the $\rho_{l}$ and $u_{g}$ scales are large and the
appearance of an unphysical $\rho_{l}u_{g}^{2}$ dynamic pressure scale could
create numerical pressure fluctuations of the same order and unphysical
pressure spikes, as the one nicely illustrated in [10] in the front part of a
suddenly accelerated small droplet. One way of avoiding this unphysical
mixture of liquid and gas quantities is to extrapolate liquid and gas pressure
and velocity in the “other” phase, as in the ghost fluid method. This
extrapolation strategy was used successfully in [11].

It may be argued that a way to avoid this numerical diffusion of liquid and
gas quantities is to advect the volume fraction and the conserved quantities
that depend on it (density, momentum and energy) in a consistent manner. In
incompressible flow in which we specialise in this paper, it means that the
volume fraction and momentum or velocity must be advected in the same way.
This is equivalent to request that the discontinuity of the Heaviside function
$H$, marking the phase transition, should be advected at the same speed as the
discontinuity in momentum. This can be expressed by the following consistency
requirement: if momentum is initially exactly proportional to volume fraction,
it should remain so after advection. We call such a method VOF-consistent.

To satisfy this requirement the idea is to solve the advection equation for
momentum with the same numerical scheme that is used for the VOF color
function. This consistency property minimizes the unphysical transfer of
momentum from one phase to another due to the differences in the numerical
schemes used. The consistency is especially important when dealing with fluids
with a large density contrast where a small numerical momentum transfer from
the dense phase to the light phase results in large numerical errors in the
velocity field which in turn creates numerical instabilities.

In this work, we present a modification of the classical momentum-preserving
scheme proposed by [1] for the case of a staggered grid and VOF method. We
also focus more on the consistency requirement between VOF and momentum
advection, and only partially conserve momentum in the scheme. This will
illustrate mostly the benefit of using a consistent advection scheme.

The paper is organized as follow: the second section deals with the continuum
mechanics formulation for incompressible flow and sharp interfaces. Section 3
describes our numerical method, starting with an overview of already-known
methods for spatial discretization, time-stepping, and VOF advection. We
continue with the new momentum advection-VOF-consistent method. Section 4 is
devoted to tests of the method, followed by a conclusion.

Among the authors, Gretar Tryggvason, Ruben Scardovelli, Yue Stanley Ling and
Stephane Zaleski have been involved in the construction of the base of the
ParisSimulator VOF and Front-Tracking code that was used to implement and test
the ideas in this paper. ParisSimulator is itself based on a Front-Tracking
code developed by Gretar Tryggvason, Sadegh Dabiri and Jiacai Lu. Daniel
Fuster was involved in the development of the Momentum Conserving method with
consistent VOF advection, with help from Tomas Arrufat, Leon Malan and Yue
Stanley Ling. The falling raindrop testing and the corresponding figures were
done by Tomas Arrufat for the CIAM method and Sagar Pal for the WY method. The
atomisation testing was done by Marco Crialesi-Esposito.

##  2 Navier–Stokes equations with interfaces

We model flows with sharp interfaces defined implicitly by a characteristic
function $H({\bf x},t)$ defined such that fluid 1 corresponds to $H=1$ and
fluid 2 to $H=0$. The viscosities $\mu$ and densities $\rho$ are calculated as
an average

$\mu=\mu_{1}H+\mu_{2}(1-H)\,,\qquad\rho=\rho_{1}H+\rho_{2}(1-H)\,.$ (1)

There is no phase change so the interface, almost always a smooth
differentiable surface $S$, advances at the speed of the flow, that is
$V_{S}={\bf u}\cdot{\bf n}$ where ${\bf u}$ is the local fluid velocity and
${\bf n}$ a unit normal vector perpendicular to the interface. Equivalently
the interface motion can be expressed in weak form

$\partial_{t}H+{\bf u}\cdot{\nabla}H=0\,,$ (2)

which expresses the fact that the singularity of $H$, located on $S$, moves at
velocity $V_{S}={\bf u}\cdot{\bf n}$. For incompressible flows, which we will
consider in what follows, we have

$\nabla\cdot{\bf u}=0\,.$ (3)

The Navier–Stokes equations for incompressible, Newtonian flow with surface
tension may conveniently be written in operator form

$\partial_{t}(\rho{\bf u})={\cal L}(\rho,{\bf u})-{\nabla}p$ (4)

where ${\cal L}={\cal L}_{\rm conv}+{\cal L}_{\rm diff}+{\cal L}_{\rm
cap}+{\cal L}_{\rm ext}$ so that the operator ${\cal L}$ is the sum of
advective, diffusive, capillary force and external force terms. The first two
terms are

${\cal L}_{\rm conv}=-{\nabla}\cdot(\rho{\bf u}{\bf u})\,,\qquad{\cal L}_{\rm
diff}={\nabla}\cdot{\bf D}\,,$ (5)

where ${\bf D}$ is a stress tensor whose expression for incompressible flow is

${\bf D}=\mu\left[\nabla{\bf u}+(\nabla{\bf u})^{T}\right]\,,$ (6)

where $\mu$ is computed from $H$ using (1). The capillary term is

${\cal L}_{\rm cap}={\sigma\kappa\delta_{S}{\bf
n}}\,,\qquad\kappa=1/R_{1}+1/R_{2}\,,$ (7)

where $\sigma$ is the surface tension, ${\bf n}$ is the unit normal
perpendicular to the interface, $\kappa$ is the sum of the principal
curvatures and $\delta_{S}$ is a Dirac distribution concentrated on the
interface. We assume a constant surface tension value $\sigma$. Finally ${\cal
L}_{\rm ext}$ represents external forces such as gravity.

##  3 Method

###  3.1 Spatial discretization

Figure 1: Representation of the staggered spatial discretisation. The pressure
$p$ is assumed to be known at the center of the control volume outlined by a
thick solid line. The horizontal velocity component $u_{1}=u$ is stored in the
middle of the left and right edges of this control volume and the vertical
velocity component $u_{2}=v$ in the middle of the top and bottom edges.

We assume a regular cubic or square grid. This can be easily generalized to
rectangular or cuboid grids, and with some efforts to quadtree and octree
grids. We also use staggered velocity and pressure grids. This makes our
method more complex than it would be on a collocated grid.

The staggered grid is represented in Figure 1. In a staggered grid
discretization of the advection equation the control volumes of the velocity
components $u_{1}$ and $u_{2}$ are shifted with respect to the control volume
surrounding the pressure $p$. The use of staggered control volumes has the
advantage of suppressing neutral modes often observed in collocated methods
but leads to more complex discretizations (see [12] for a more detailed
discussion.) The control volumes for $u_{1}$ and $u_{2}$ are shown on Figure
2. This type of staggered representation is easily generalized to three
dimensions. In what follows we shall use these control volumes for the
velocity or momentum components.

Using the staggered grid leads to a compact expression for the continuity
equation (3)

${u_{1;i+1/2,j,k}-u_{1;i-1/2,j,k}\over\Delta
x}+{u_{2;i,j+1/2,k}-u_{2;i,j-1/2,k}\over\Delta
y}+{u_{3;i,j,k+1/2}-u_{3;i,j,k-1/2}\over\Delta z}=0,$ (8)

(a) (b)

Figure 2: The control volumes for the $u_{1}=u$ and the $u_{2}=v$ velocity
components are displaced half a grid cell to the right (horizontal velocities)
and to the top (vertical velocities). Here the indices show the location of
the stored quantities. Thus, half indices indicate those variables stored at
the edges of the pressure control volumes.

In what follows, we shall use the notation $f=m\pm$, with the integer index
$m=1,2,3$, to note the face of any control volume located in the positive or
negative Cartesian direction $m$, and ${\bf n}_{f}$ for the normal vector of
face $f$ pointing outwards of the control volume. On a cubic grid the spatial
step is $\Delta x=\Delta y=\Delta z=h$ so the continuity equation becomes

${\nabla}^{h}\cdot{\bf u}=\sum_{m=1}^{3}(u_{m+}+u_{m-})/h=0\,,$ (9)

where $u_{f}=u_{m\pm}={\bf u}\cdot{\bf n}_{f}$ is the velocity normal to face
$f$. The discretization of the interface location is performed using a VOF
method. VOF methods typically attempt to solve approximately equation (2)
which involves the Heaviside function $H$, whose integral in the cell $\Omega$
indexed by $i,j,k$ defines the volume fraction $C_{i,j,k}$ from the relation

$h^{3}\,C_{i,j,k}=\int_{\Omega}H{\rm d}{\bf x}\,.$ (10)

$C_{i,j,k}$ represents the fraction of the cell labelled by $i,j,k$ filled
with fluid 1, taken to be the reference fluid. It is worth noting that in the
staggered grid setup, the control volume for $p$ is also the control volume
for other scalar quantities, such as $\rho$ and $H$.

###  3.2 Time Marching

The volume fraction field is updated as

$C^{n+1}=C^{n}+{\cal L}_{\rm VOF}(C^{n},{\bf u}^{n}\tau/h)\,,$ (11)

where ${\cal L}_{VOF}$ represents the operator that updates the Volume of
Fluid data given the velocity field. Once volume fraction is updated, the
velocity field is updated in a couple of steps. A projection method is first
used, in which a provisional velocity field ${\bf u}^{*}$ is computed

$\rho^{n+1}{\bf u}^{*}=\rho^{n}{\bf u}^{n}+\tau{\cal L}^{h}_{\rm
conv}(\rho^{n},{\bf u}^{n})+\tau\left[{\cal L}^{h}_{\rm diff}(\mu^{n},{\bf
u}^{n})+{\cal L}^{h}_{\rm cap}(C^{n+1})+{\cal L}^{h}_{\rm
ext}(C^{n+1})\right]$ (12)

It goes without saying that the above operators depend on the discretization
steps $\tau$ and $h$ as well as the fluid parameters. The discussion of the
${\cal L}^{h}_{\rm conv}$ operator is the main point of this paper. In the
second step, the projection step, the pressure gradient force is added to
yield the velocity at the new time step

${\bf u}^{n+1}={\bf u}^{*}-\frac{\tau}{\rho^{n+1}}\nabla^{h}p\,.$ (13)

The pressure is determined by the requirement that the velocity at the end of
the time step must be divergence free

${\nabla}^{h}\cdot{\bf u}^{n+1}=0\,,$ (14)

which leads to a Poisson-like equation for the pressure

${\nabla}^{h}\cdot\frac{\tau}{\rho^{n+1}}\nabla^{h}p={\nabla}^{h}\cdot{\bf
u}^{*}\,.$ (15)

###  3.3 Volume-Of-Fluid

We follow the usual VOF method and will only detail the necessary steps to
illustrate the momentum advection method based on the VOF method in the next
section. Moreover we will use in this section a rescaling of the space and
time variables so that the cell size is 1, and the time step is also 1. All
velocities are then rescaled to $u^{\prime}=u\tau/h$. Because of this space
rescaling and in these new units, $C_{i,j,k}$ is also the measure of the
volume of reference fluid in cell $i,j,k$.

####  3.3.1 Volume initialisation

Before any VOF interface tracking is performed, the field of $C_{i,j,k}$
values must be initialized. We use the Vofi library described in [13] and
[14]. This allows a high-accuracy numerical integration of the measure of the
fluid volumes.

####  3.3.2 Normal vector determination

The VOF method proceeds in two steps, reconstruction and advection. In the
reconstruction step, one attempts to obtain the “best” possible representation
of the interface using Volume of Fluid data $C_{i,j,k}$. The best
representation is usually the one with the smallest reconstruction error for
objects such as circles or ellipses, but issues of continuity and robustness
at low resolution may also be considered. In the work reported here we use a
standard method described in [12]. In each cell cut by the interface one first
determines the interface normal vector ${\bf n}$, then one solves the problem
of finding a plane perpendicular to ${\bf n}$ under which one finds exactly
the volume $C_{i,j,k}$. Normal vector determination is performed using the MYC
method described in [12].

####  3.3.3 Plane constant determination

Once the interface normal vector ${\bf n}$ is determined, a new, colinear
normal vector noted ${\bf m}$ and having unit $L_{1}$ norm is deduced from
${\bf n}$, so that $|m_{x}|+|m_{y}|+|m_{z}|=1$. Considering the volume
$V=C_{i,j,k}$ in cell $i,j,k$ the plane constant $\alpha$ is defined so that
the plane

${\bf m}\cdot{\bf x}=\alpha$ (16)

cuts exactly a volume $V$ in the cell. Typically $\alpha$ is determined by the
resolution of a cubic equation as described in [15].

####  3.3.4 General split-direction advection

Once the reconstruction has been performed at time $t_{n}$, it is used to
obtain the approximate position of the interface, and the volumes $C_{i,j,k}$
at time $t_{n+1}$. The following discussion of momentum advection is based on
two VOF advection methods, (Figure 4) Lagrangian Explicit and Weymouth and Yue
advection. We first describe the common features of these two methods.

After addition and subtraction of a term proportional to the velocity
divergence, equation (2) leads to

$\partial_{t}H+\nabla\cdot({\bf u}H)=H\,\nabla\cdot{\bf u}\,.$ (17)

A discrete representation of this equation is obtained by integrating the
advection equation in time and space

${C_{i,j,k}^{n+1}-C_{i,j,k}^{n}}=-\sum_{\rm{faces}\,f}F^{(c)}_{f}+\int_{t_{n}}^{t_{n+1}}{\rm
d}t\int_{\Omega}H\,\nabla\cdot{\bf u}\,{\rm d}{\bf x}\,,$ (18)

where the first term on the right-hand side is the sum over faces $f$ of cell
$i,j,k$ of the fluxes $F^{(c)}_{f}$ of ${\bf u}H$. Obviously the “compression”
term on the right-hand side disappears for incompressible flow, however it is
useful to keep it in view of its utility in split-advection methods. In the
previous equation $F_{f}^{(c)}$ is defined by

$F_{f}^{(c)}=\int_{t_{n}}^{t_{n+1}}{\rm d}t\int_{f}u_{f}({\bf x},t)H({\bf
x},t)\,{\rm d}{\bf x}\,,$ (19)

where $u_{f}={\bf u}\cdot{\bf n}_{f}$ and ${\bf n}_{f}$ is the unit normal
vector of face $f$, pointing outwards of cell $i,j,k$.

Once an approximation for the evolution of ${\bf u}({\bf x},t)H({\bf x},t)$
during the time step is chosen, a four-dimensional integral remains to be
computed in equation (19). In what follows we describe the approximations
used. However, all the methods we describe here are directionally split. The
methods are also designed to preserve the property that $0\leq C_{i,j,k}\leq
1$ which we call $C$-bracketing. It is important to preserve $C$-bracketing
despite the various approximations used, in order to avoid the arbitrary
addition or removal of mass.

Directional splitting results in the breakdown of equation (18) into three
equations

${C_{i,j,k}^{n,m+1}-C_{i,j,k}^{n,m}}=-F^{(c)}_{m-}-F^{(c)}_{m+}+c_{m}\partial_{m}^{h}u_{m}\,,$
(20)

where $m=1,2,3$, $C_{i,j,k}^{n,1}=C_{i,j,k}^{n}$,
$C_{i,j,k}^{n,4}=C_{i,j,k}^{n+1}$, the face labelled “$m-$” is the “left” face
in direction $m$ with $F^{(c)}_{m-}\geq 0$ if the flow is locally from “right”
to “left”, and finally the face “$m+$” is the “right” face in direction $m$
with $F^{(c)}_{m+}\geq 0$ if the flow is from “left” to “right”. We have also
approximated the compression term in (18) by

$\int_{t_{n}}^{t_{n+1}}{\rm d}t\int_{\Omega}H\partial_{m}u_{m}{\rm d}{\bf
x}\simeq c_{m}\partial_{m}^{h}u_{m}\,,$ (21)

with no implicit summation rule. In the RHS of (20) and 21) the flux terms
$F_{f}^{(c)}$ and the partial derivative $\partial_{m}u_{m}$ must be evaluated
with the same discretized velocities. In particular, $\partial_{m}^{h}u_{m}$
is a finite difference or finite volume approximation of the spatial
derivative of the $m$th component of the velocity vector in direction $m$, and
the “compression coefficient” $c_{m}$ approximates the color fraction. Its
exact expression is dependent on the version of the directionally split
advection method used and will be described below. The value of $c_{m}$ is
chosen so as to preserve the $C$-bracketing condition. Note that because
$c_{m}$ may not be the same coefficient along the three Cartesian directions,
even if the flow is incompressible the sum
$\sum_{m}c_{m}\partial_{m}^{h}u_{m}$ of the split compression terms is not
necessarily vanishing.

Each application of equation (20) defines an advection substep. After each
substep, a new reconstruction of the interface is performed with the updated
volumes $C_{i,j,k}^{n,m+1}$, computing the normal ${\bf m}$ and the $\alpha$
constant. This reconstruction is used to estimate the fluxes $F^{(c)}_{f}$ in
the next split advection substep (20).

####  3.3.5 Lagrangian Explicit advection

The Lagrangian Explicit / Calcul d’Interface Affine par Morceaux CIAM
advection method is described in what follows. “Calcul d’Interface Affine par
Morceaux” is French for “Piecewise Linear Interface Calculation” (PLIC) but
while PLIC refers to a generic VOF method with a piecewise linear
reconstruction step, CIAM refers to a specific type of advection method first
described in the archival literature in [16] and classified as the “LE” method
in [17]. CIAM advection is most naturally explained as a Lagrangian transport
of the Heaviside function. In each cell, the reconstruction step has defined a
planar interface. Lagrangian transport consists in moving each point of this
interface by a velocity field ${\bf u}^{\prime}$ defined as follows: 1) in
each sub time step $m$ of split advection the velocity field is in direction
$m$, and depends only on the component $x_{m}$ of the space variable; 2) the
velocity field is interpolated between the face velocities $u_{m-}$ on the
left face and $u_{m+}$ on the right face so that
$u_{m}^{\prime}(x_{m})=-u_{m-}(1-x_{m})+u_{m+}x_{m}$ where the origin of the
coordinate $x_{m}$ was chosen on the left face. (As said before, space has
been rescaled so the cell size is unity and recall $u_{f}={\bf u}\cdot{\bf
n}_{f}$.) The exact transport by this velocity field is
$dx_{m}/dt=u_{m}^{\prime}(x_{m})$ but we use a first-order, explicit
approximation as $x_{m}^{n+1}-x_{m}^{n}=u_{m}^{\prime}(x_{m}^{n})$. (Time was
rescaled so the time step is unity and disappears from the equations.) Using
the approximation for the interpolated velocity we get

$x_{m}^{n+1}=-u_{m-}+(1+u_{m+}+u_{m-})\,x^{n}_{m}\,.$ (22)

This relation gives the advected points of the interface as a function of the
original points. It is used to advect the faces and the planar interfaces as
shown on Figure 3, for a two-dimensional case. After advection, at most three
interface fragments are found in each cell. The resulting volumes (areas in
2D) are shown in Figure 4 and are labelled $V_{1}$, $V_{2}$ and $V_{3}$.
Eventually

Figure 3: Advection of the AB interface segment. Face velocities at the left
face, $u_{i-1/2,j}=-u_{1-}$, and the right face, $u_{i+1/2,j}=u_{1+}$, are
used to create an interpolated linear velocity field, so that point $A$ is
advected to point $A^{\prime}$ at the velocity $u_{i-1/2,j}$, point $B$ to
$B^{\prime}$ at the velocity $u_{i+1/2,j}$, and a point along AB at a velocity
linearly interpolated between $u_{i-1/2,j}$ and $u_{i+1/2,j}$. Figure 4:
Formation of the volumes $V_{1},V_{2}$ and $V_{3}$ by Lagrangian advection in
the $x$ direction. Left: initial reconstruction with the horizontal velocities
on the faces of the central cell. Right: segments and volumes $V_{i}$ (see
text) after Lagrangian advection (for simplicity it has been assumed
$u_{i-3/2,j}=u_{i+3/2,j}=0$)

the color function is given by the sum of three contributions

$C_{i,j,k}^{n,m+1}=V_{1}+V_{2}+V_{3}\,.$ (23)

It is seen that the transform (22) compresses distances by a factor
$1+\partial^{h}_{m}u_{m}=1+u_{m+}+u_{m-}$. Thus the measure of any region of
space is scaled by this amount under advection. The interface is then
transfomed in three steps: 1) compress all cells by the factor above; 2) move
the cells by $u_{m-}$; 3) the moved cells now overlap the original cells. The
volume corresponding to the overlap of the cell with itself is $V_{2}$, while
$V_{1}$ and $V_{3}$ correspond to the overlap with the moved cells from the
left and right, respectively. The geometrical interpretation of the Lagrangian
Explicit advection of Figs. 3 and 4 and the definition (19) of $F^{(c)}_{f}$
leads to the following correspondence between the fluxes and the volume
contributions $V_{i}$. To start with an example, for the central cell of Fig.
4 the flux on the left face is from left to right, since
$u_{1-}=-u_{i-1/2,j}<0$. This means that
$V_{1}=-F^{(c)}_{1-,i}=F^{(c)}_{1+,i-1}>0$ while $V_{3}$ of the left cell,
that is $V_{3,i-1}$, is zero corresponding to an absence in that cell. In all
cases the fluxes accross a face are equal to the $V_{1}$ or $V_{3}$ of one
cell or the other adjacent to a face. The final expression of the split step
is then

$C_{i,j,k}^{n,m+1}=C_{i,j,k}^{n,m}(1+u_{m+}+u_{m-})-F^{(c)}_{m-}-F^{(c)}_{m+}\,,$
(24)

which shows that the constant $c_{m}$ in the compression term is
$C_{i,j,k}^{n,m}$ while the approximation of the derivative is
$\partial_{m}^{h}u_{m}=u_{m+}+u_{m-}$ .

It is interesting to note that using three Lagrangian advections in a sequence
does not result in volume conservation to machine accuracy. Indeed the
summation of the three substeps (20) results in

${C_{i,j,k}^{n+1}-C_{i,j,k}^{n}}=-\sum_{\rm{faces}\,f}F^{(c)}_{f}+\sum_{m=1}^{3}C_{i,j,k}^{n,m}(u_{m+}+u_{m-})\,.$
(25)

While the flux terms cancel upon integration over the domain, the sum of the
compressive terms does not vanish since $C_{i,j,k}^{n,m}$ changes during the
three substeps. We recover eq. (20) with

$c_{m}=C_{i,j,k}^{n,m}$ (26)

Figure 5: Eulerian flux representation for advection in the $x$ direction.
Left: same initial reconstruction of Fig. 4 with the horizontal velocities on
the faces of the central cell. Right: fluxes, or volumes $V_{1}$ and $V_{3}$,
are calculated directly from the interface reconstruction in each cell

####  3.3.6 Weymouth and Yue advection

Another method that we shall describe is the exactly mass-conserving method of
Weymouth and Yue (WY) described in ref. [18]. In that method the coefficient
of the compression term (21) is independent of the direction $m$ so that
$c_{m}=c$. The coefficient $c$ is defined as

$c=\Theta(C_{i,j,k}^{n}-1/2)$ (27)

where $\Theta$ is the one-dimensional Heaviside function (defined on the real
axis). That is $c=0$ if $C_{i,j,k}^{n}<1/2$ and $c=1$ if $C_{i,j,k}^{n}\geq
1/2$. The fluxes $F^{(c)}_{f}$ are also defined differently. The left fluxed
volume in direction $m=1$ is equal to the volume fraction $F_{1-}^{(c)}$ in a
cuboid of width $u_{i-1/2,j,k}$ (recall that $\tau=1$) adjacent to the left
face $f=1-$. This fluxed volume corresponds to “Eulerian Implicit” (EI)
advection in the terminology of [17] and is represented by the volume $V_{1}$
on Figure 5. Using these definitions, Weymouth and Yue were able to show in
ref. [18] that the final result obeys $C$-bracketing (see Section 3.3.4).

Using this advection method results in volume conservation at machine
accuracy. Indeed the summation of the three substeps (20) results in

${C_{i,j,k}^{n+1}-C_{i,j,k}^{n}}=-\sum_{\rm{faces}\,f}F^{(c)}_{f}+c\sum_{m=1}^{3}\partial_{m}^{h}u_{m}$
(28)

Since $\sum_{m=1}^{3}\partial_{m}^{h}u_{m}=\sum_{m=1}^{3}(u_{m+}+u_{m-})$ is
the simple finite-volume expression for $\nabla\cdot{\bf u}$, it disappears
and mass is conserved at the accuracy with which condition (3) is satisfied.

####  3.3.7 Clipping

The algorithm that has been coded involves a number of additional steps
designed to avoid unwanted effects of arithmetic floating point round-off
error. The most important one is clipping: at the end of each directional
advection, the values of $C_{i,j,k}$ are reset so that $C_{ijk}$ is set to $0$
is $C_{ijk}<\epsilon_{c}$ and $C_{ijk}$ is set to $1$ is
$C_{ijk}>1-\epsilon_{c}$. When there is no surface tension the choice
$\epsilon_{c}=10^{-12}$ works well. Otherwise $\epsilon_{c}=10^{-8}$ gives
more stable results with smoother interface shapes. This stronger clipping is
a necessity for some simulations with WY, while the CIAM method can survive
$\epsilon_{c}=10^{-12}$ with the diffusive interpolation schemes described
below. As a matter of fact we observe that WY produces many more “wisps”, i.e.
cells with tiny values of $1-C_{i,j,k}$ inside the liquid or $C_{i,j,k}$
inside the gas. We have not yet been able to determine the origin of this need
for a more forceful clipping with WY, but it could be related to the fact that
the CIAM method has a geometrical interpretation, while WY is intrisically
algebraic in nature.

###  3.4 Momentum-advection methods

####  3.4.1 Advection of a generic conserved quantity

Consider the advection of a generic conserved quantity $\phi$ by a continuous
velocity field

$\partial_{t}\phi+\nabla\cdot(\phi{\bf u})=0.$ (29)

We assume that $\phi$ is smoothly varying except on the interface where it may
be discontinuous. Indeed finding a correct scheme for the advection of this
discontinuity, at the same speed as the advection of the volume fraction, is
the goal of the present study. The smoothness of the advected quantity away
from the interface is verified for the density $\rho$, the momentum $\rho{\bf
u}$ or the internal energy $\rho e$. Integrating the advection equation in
time one obtains

${\phi_{i,j,k}^{n+1}-\phi_{i,j,k}^{n}}=-\sum_{\rm{faces}\,f}F^{(\phi)}_{f},$
(30)

where the sum of the right-hand side is the sum over faces $f$ of cell $i,j,k$
of the fluxes $F^{(\phi)}_{f}$ of $\phi$. In equations, $F_{f}^{(\phi)}$ is
defined in a similar manner to equation (19) for the color function fluxes.

$F_{f}^{(\phi)}=\int_{t_{n}}^{t_{n+1}}{\rm d}t\int_{f}u_{f}({\bf
x},t)\phi({\bf x},t){\rm d}{\bf x},$ (31)

where

$u_{f}={\bf u}\cdot{\bf n}_{f},$ (32)

and ${\bf n}_{f}$ is the unit normal vector of face $f$, pointing outwards of
cell $i,j,k$. In order to “extract” the discontinuity the flux can be
expressed as

$F_{f}^{(\phi)}=\int_{t_{n}}^{t_{n+1}}{\rm
d}t\int_{f}[u_{f}H\phi+u_{f}(1-H)\phi]{\rm d}{\bf x},$ (33)

which can be re-expressed in terms of the fluxes $F^{(c)}$ in eq. (19) as

$F_{f}^{(\phi)}=\bar{\phi}_{1}\int_{t_{n}}^{t_{n+1}}{\rm d}t\int_{f}u_{f}H{\rm
d}{\bf x}+\bar{\phi}_{2}\int_{t_{n}}^{t_{n+1}}{\rm d}t\int_{f}u_{f}(1-H){\rm
d}{\bf x}$ (34)

where

$\bar{\phi}_{s}=\frac{\int_{t_{n}}^{t_{n+1}}{\rm
d}t\int_{f}u_{f}H_{s}\phi\,{\rm d}{\bf x}}{\int_{t_{n}}^{t_{n+1}}{\rm
d}t\int_{f}u_{f}H_{s}\,{\rm d}{\bf x}},$ (35)

and $H_{1}=H$, $H_{2}=1-H$. This can be expressed as

$F_{f}^{(\phi)}=\bar{\phi}_{1}F_{f}^{(c)}+\bar{\phi}_{2}F_{f}^{(1-c)},$ (36)

where $F_{f}^{(c)}$ is the flux of $H$ defined in (19) and $F_{f}^{(1-c)}$ is
the flux obtained by replacing $H$ by $1-H$ in (19).

####  3.4.2 Cloning the tracers

When a cell is cut by the interface, and the field $\phi$ is non-smooth, it
becomes difficult to estimate the integrals in (34). A possibility is to
define two new fields $\phi_{1,2}$ so that $\phi_{s}=\phi$ inside phase $s$.
Then

$\phi=H\phi_{1}+(1-H)\phi_{2}.$ (37)

In the discretized equation, two discrete fields $\phi_{s,i,j,k}$ are defined
over the entire domain $\Omega$. This is more costly in memory usage but
simplifies considerably the computation of the averages in Equation (34). The
two equations (29) and (2) are now replaced by three equations, the same
volume fraction equation (2) and the two equations

$\displaystyle\partial_{t}\phi_{1}+\nabla\cdot(\phi_{1}{\bf u})$
$\displaystyle=0$ , (38)
$\displaystyle\partial_{t}\phi_{2}+\nabla\cdot(\phi_{2}{\bf u})$
$\displaystyle=0$ . (39)

The three equations (2,38-39) now imply (29). This addition of a pair of
“cloned” variables to help deal with large variations of density is similar to
the methods used for the resolution of the momentum and energy equations for
compressible flow. For example Saurel and Abgrall used two density, momentum
and energy variables in ref. [7], leading to their seven-equations model,
while Allaire, Clerc and Kokh use two density variables in [8], leading to
their five-equations model. The addition of a cloned tracer variable in
incompressible isothermal flow was also implemented by Popinet in the
“Basilisk” code [19].

####  3.4.3 Advection of the density field

Let us now apply the above considerations to the advection of the density
field by a general continuous velocity field. Although in the incompressible
case the density follows trivially the color function, it will allow us to
introduce an important point about the equivalence of tracer and VOF
advection. The density $\rho({\bf x},t)$ obeys eq. (29) with $\phi=\rho$. The
advection of density may be made consistent with the advection of the color
function by using fluxes $F_{f}^{(c)}$ and estimating the average occuring in
equation (35), $\bar{\phi}_{1}=\bar{\rho}$ from averages over control volumes
or approximations on faces.

The case relevant here is that of an incompressible velocity field, with
uniform density in each phase. Then $\rho$ is exactly proportional to $H$ (Eq.
1). The constancy of $\rho$ allows to extract it trivially from the integrals
in eq. (34) and to obtain exactly $\bar{\rho}_{s}=\rho_{s}$. Then the flux of
$\rho$ is

$F_{f}^{(\rho)}=\rho_{1}F_{f}^{(c)}+\rho_{2}F_{f}^{(1-c)}.$ (40)

Using this flux definition for $\rho$, and any VOF method for the fluxes of
the color function, one obtains a conservative method for $\rho$ since eq.
(30) evolves $\rho$ as a difference of fluxes. That is, one conserves the
total mass and
$\sum_{i,j,k}\rho^{(n)}_{i,j,k}=\sum_{i,j,k}\rho^{(n+1)}_{i,j,k}$. However
this result is not consistent with the advection of the color function in the
CIAM case, since as we have seen above, CIAM does not conserve volumes
exactly, that is $\sum_{i,j,k}C^{(n)}_{i,j,k}=\sum_{i,j,k}C^{(n+1)}_{i,j,k}$
is not in general true. As a result the advection of $\rho$ is not consistent
with the advection of $C$.

The paradox may be resolved if one notices that the compression term is
missing in the equation for the advection of $\rho$. One should keep the
compression term, which is not changing the equation in the case of
incompressible velocity fields, so that the equation for the conserved
quantity $\phi$ or $\rho$ becomes

$\partial_{t}\phi+\nabla\cdot(\phi{\bf u})=\phi\,\nabla\cdot{\bf u}$ (41)

It is then possible to define the evolution of $\rho=\phi$ through a sequence
of directionally split operations exactly equivalent to the operations
performed on the color function.

${\phi_{i,j,k}^{n,m+1}-\phi_{i,j,k}^{n,m}}=-F^{(\phi)}_{m-}-F^{(\phi)}_{m+}+(\tilde{\phi}_{1}^{m}c^{(1)}_{m}+\tilde{\phi}_{2}^{m}c^{(2)}_{m})\partial_{m}^{h}u_{m}$
(42)

where $F^{(\phi)}_{m\pm}$ is defined as above, and

$\tilde{\phi}_{s}^{m}=\frac{\int_{t_{n}}^{t_{n+1}}{\rm d}t\int_{\Omega}\phi
H_{s}\partial_{m}^{h}u_{m}\,{\rm d}{\bf x}}{\int_{t_{n}}^{t_{n+1}}{\rm
d}t\int_{\Omega}H_{s}\partial_{m}^{h}u_{m}\,{\rm d}{\bf x}}$ (43)

and where $c^{(1)}_{m}=c_{m}$ is the compression coefficient computed as in
the VOF advection from $C_{i,j,k}^{n,m}$ and $c^{(2)}_{m}=1-c_{m}$ corresponds
to the symmetric color fraction $1-C_{i,j,k}^{n,m}$. Specifically for $\rho$
this gives

${\rho_{i,j,k}^{n,m+1}-\rho_{i,j,k}^{n,m}}=-F^{(\rho)}_{m-}-F^{(\rho)}_{m+}+C_{m}^{(\rho)},$
(44)

where the fluxes are given by (40) and the central compression part is

$C_{m}^{(\rho)}=(\rho_{1}c^{(1)}_{m}+\rho_{2}c^{(2)}_{m})\partial_{m}^{h}u_{m}$
(45)

with no implicit summation on $m$. The $c_{m}^{(s)}$ values are given by (26)
or (27). It is interesting to note that for the WY method, the compression
terms eventually cancel and mass is conserved at the same accuracy as the
discrete incompressibility condition $\partial_{m}^{h}u_{m}=0$ is verified.

####  3.4.4 Momentum advection: basic expressions

Momentum advection can be performed following the general outline above. Each
momentum component can be treated by considering the transport of the
conserved quantity $\phi=\rho u_{q}$ where $q$ is the component index.
Reworking equations (33,34) and (35), with this definition, we obtain for the
weighted averages $\bar{\phi}_{s}=\overline{\rho u_{q}}_{s}$ in phase $s$ the
expression

$\overline{\rho u_{q}}_{s}=\rho_{s}\bar{u}_{q,s}$ (46)

where

$\bar{u}_{q,s}=\frac{\int_{t_{n}}^{t_{n+1}}{\rm
d}t\int_{f}u_{q}u_{f}H_{s}\,{\rm d}{\bf x}}{\int_{t_{n}}^{t_{n+1}}{\rm
d}t\int_{f}u_{f}H_{s}\,{\rm d}{\bf x}}$ (47)

We term $\bar{u}_{q,s}$ the “advected interpolated velocity” and we will
explain below how these “advected” velocities are computed. Thus the evolution
of the momentum is given by

$\displaystyle{(\rho_{i,j,k}u_{q;i,j,k})^{n,m+1}-(\rho_{i,j,k}u_{q;i,j,k})^{n,m}}$
$\displaystyle=$ $\displaystyle-F^{(\rho u)}_{m-}-F^{(\rho
u)}_{m+}+(\rho_{1}\tilde{u}_{q,1}^{m}c^{(1)}_{m}$ $\displaystyle+$
$\displaystyle\rho_{2}\tilde{u}_{q,2}^{m}c^{(2)}_{m})\partial_{m}^{h}u_{m}$
(48)

where

$F^{(\rho
u)}_{f}=\rho_{1}\bar{u}_{q,1}F^{(c)}_{f}+\rho_{2}\bar{u}_{q,2}F^{(1-c)}_{f},$
(49)

and the “central interpolated velocity” corresponding to the
$\tilde{\phi}_{s}$ are

$\tilde{u}_{q,s}^{m}=\frac{\int_{t_{n}}^{t_{n+1}}{\rm
d}t\int_{\Omega}u_{q}^{n,p}H_{s}\partial_{m}^{h}u_{m}\,{\rm d}{\bf
x}}{\int_{t_{n}}^{t_{n+1}}{\rm
d}t\int_{\Omega}H_{s}\partial_{m}^{h}u_{m}\,{\rm d}{\bf x}}$ (50)

where for CIAM $p=m$ and for WY $p=1$. Thus
$\tilde{u}_{q,s}^{m}=\tilde{u}_{q,s}$ remains constant while the various steps
$m=1,..,3$ of directional splitting are performed. Note that we otherwise omit
the superscript $m$ for $\bar{u}_{q}$ and $\tilde{u}_{q}$ to avoid too complex
notations. Notice that “cloning” the advected velocities $\bar{u}_{q,1}$ and
$\bar{u}_{q,2}$ would make it easier to advect a velocity field with a jump on
the interface. However in viscous flow without phase change the velocity is
continuous on the interface, and to avoid an excessively complicated method we
approximate the velocity field as continuous and we choose approximations of
the “advected interpolated velocity” such that
$\bar{u}_{q}=\bar{u}_{q,1}=\bar{u}_{q,2}$ and “central interpolated velocity”
such that $\tilde{u}_{q}=\tilde{u}_{q,1}=\tilde{u}_{q,2}$. An important
simplification is then

$F^{(\rho u)}_{f}=\bar{u}_{q}F^{(\rho)}_{f}$ (51)

(which is the central equation in this development) and thus

${(\rho_{i,j,k}u_{q;i,j,k})^{n,m+1}-(\rho_{i,j,k}u_{q;i,j,k})^{n,m}}=-\bar{u}_{q}F^{(\rho)}_{m-}-\bar{u}_{q}F^{(\rho)}_{m+}+\tilde{u}_{q}^{m}C_{m}^{(\rho)},$
(52)

where the density fluxes are defined in (40) and the central term $C^{(\rho)}$
is defined in (45). In the above expression the weighted average velocities
$\bar{u}_{q}$ are computed using (47) on the corresponding left face $m-$ or
right face $m+$.

The scheme above must be combined with an advection scheme for the color
fraction. When the CIAM scheme is used, the $C_{m}^{(\rho)}$ term in (52) must
be kept to ensure consistency with VOF advection. This term is computed using
the above definition (26). This term does not cancel when the final momentum
is computed from the three directionally-split advections and the result is
not exactly conservative. On the other hand when the WY scheme is used the
$c^{(s)}_{m}$ and $\tilde{u}_{q}^{m}$ factors are independent of $m$ and

$C_{m}^{(\rho)}=[\rho_{1}c+\rho_{2}(1-c)]\partial_{m}^{h}u_{m}$ (53)

so provided the velocity field is incompressible (that is
$\sum_{m=1}^{3}\partial_{m}^{h}u_{m}=0$) summing the three split advection
expressions (52) one obtains a cancellation of the central terms and

${(\rho_{i,j,k}u_{q;i,j,k})^{n,4}-(\rho_{i,j,k}u_{q;i,j,k})^{n}}=-\sum_{\rm
faces\,f}\bar{u}_{q,f}F^{(\rho)}_{f}.$ (54)

The version of our method coupled with WY advection is thus exactly
conservative. Both CIAM and WY methods can be rewritten in a slightly
different notation useful for future discussion

${\hat{g}_{q}^{n+1}-g_{q}^{n}}=-\sum_{\rm
faces\,f}\bar{u}_{q,f}F^{(\rho)}_{f}(C^{n},u_{f})+\sum_{m=1}^{3}\tilde{u}_{q}^{m}C_{m}^{(\rho)},$
(55)

where $\hat{g}_{q}^{n+1}=(\rho_{i,j,k}u_{q;i,j,k})^{n,4}$ is the $q$-component
of the final momentum, $g_{q}^{n}=(\rho_{i,j,k}u_{q;i,j,k})^{n}$ is the
initial one, and we have written explicitly the dependence of
$F^{(\rho)}_{f}(C^{n},u_{f})$ on the face velocities and on the volume
fraction field in the neighborhood of the interface.

####  3.4.5 Momentum advection: interpolations and flux limiters

The evolution of momentum as expressed in eq. (55) can be approximated either
1) away from the interface, in the bulk of the phases or 2) in the
neighborhood of the interface. In the first case the expression simplifies
considerably since the density as well as the color fraction are constant thus
$\rho_{ijk}=\rho_{0}$ and $C^{n}=C_{0}$. The spurious central term is also
removed so

${u^{*}_{q}-u_{q}^{n}}=-\sum_{\rm faces\,f}\bar{u}_{q,f}u_{f}.$ (56)

We can distinguish an “advecting” velocity $u_{f}$ and an “advected velocity
component” $\bar{u}_{q,f}$ which must be estimated in an integral relating to
face $f$. Both components are either centered on the face of the staggered set
of cells / control volumes for the component $u_{q}$ or estimated close to the
face $f$. Typical hyperbolic equation methods use interpolants to perform
these estimations. The interpolants we use here are one-dimensional and work
from the node velocities on a segment aligned with the direction of advection,
that is the direction perpendicular to face $f$. Thus the scheme in the bulk
is

${u^{*}_{q}-u_{q}^{n}}=-\sum_{\rm
faces\,f}\bar{u}_{q,f}^{\rm(advected)}u^{\rm(advecting)}_{f}.$ (57)

Near the interface we have instead

${\hat{g}_{q}^{n+1}-g_{q}^{n}}=-\sum_{\rm
faces\,f}\bar{u}_{q,f}^{\rm(advected)}F^{(\rho)}_{f}(C^{n},u_{f}^{\rm(advecting)})+\sum_{m=1}^{3}\tilde{u}_{q}^{m}C_{m}^{(\rho)},$
(58)

To estimate the advecting velocities $u_{f}^{\rm(advecting)}$ we use centered
schemes. (see Figure 1 for a visual depiction of the arrangement of the
staggered nodes.) Consider for example that the face is perpendicular to
direction $1$ so $f=1_{\pm}$. There are two cases. In the first case the
advected component is not aligned with the face normal so in the example
$q\neq 1$. Consider specifically the case $q=2$ depicted in Figure 1(b). The
staggered $u_{2}$ control volumes are centered on $i,j+1/2,k$. The face
$f=1_{\pm}$ is then centered on $i-1/2,j+1/2,k$. Then
$u_{f}^{\rm(advecting)}=u_{1_{-}}^{\rm(advecting)}=u_{1;i-1/2,j+1/2,k}$ which
is not given and has to be interpolated by

$u_{1;i-1/2,j+1/2,k}^{\rm(advecting)}=\frac{1}{2}(u_{1;i-1/2,j,k}+u_{1;i-1/2,j+1,k}).$
(59)

In the second case the advected component is aligned with the face normal so
in the example $q=1$ also depicted in Figure 1(a). The staggered $u_{1}$
control volumes are centered on $i+1/2,j,k$. The face $f=1_{\pm}$ is then
centered on $i-1,j,k$ and the interpolation is

$u_{1;i-1,j,k}^{\rm(advecting)}=\frac{1}{2}(u_{1;i-1/2,j,k}+u_{1;i+1/2,j,k}).$
(60)

Now we turn to the definition of the advected velocity values $u_{q,f}$ in
expression (58). We still take the example where the face is perpendicular to
direction $1$ so $f=1_{\pm}$ (we take the case $f=1_{-}$ to illustrate, see
Figure 6(a)) thus the advecting direction is $x_{1}$. We only need to consider
the dependence of the advected velocity on $x_{1}$ to define the
interpolation. We consider the regularly spaced one-dimensional grid of
advected velocities. In order to use lighter notations we will note
$\phi=u_{q}$ the advected velocity component, with known node values at

$\phi_{p}=u_{q}(x_{1,1_{-}}+ph)$ (61)

where $x_{1,1_{-}}$ is the coordinate of the face $f=1_{-}$. The known values
$\phi_{p}$ in the staggered grid are at nodes that are a half-integer away
from the face, which implies that the known values are

$u_{1;i-3/2,j,k}=\phi_{-3/2},\quad u_{1;i-1/2,j,k}=\phi_{-1/2},\quad
u_{1;i-1/2,j,k}=\phi_{1/2},\cdots$ (62)

when $q=1$, $x_{1,1_{-}}=ih$ and

$u_{2;i-2,j+1/2,k}=\phi_{-3/2},\quad u_{2;i-1,j+1/2,k}=\phi_{-1/2},\quad
u_{2;i,j+1/2,k}=\phi_{1/2},\cdots$ (63)

when $q=2$. For $q=2$ or 3, we have $x_{1,1_{-}}=(i-1/2)h$ so the $q=3$ case
follows easily. We need to predict $\phi_{0}$ to serve as an approximation of
$\bar{u}_{q}$ given in (47) . An interpolation function is now defined that
predicts this value as a function of the four nearest points, and in an upwind
manner based on the sign of the advecting velocity $u_{f}$, given by centered
interpolations. The value at the face is given by

$\phi_{0}=f(\phi_{-3/2},\phi_{-1/2},\phi_{1/2},\phi_{3/2},{\rm sign}(u_{f}))$
(64)

In this study we have extensively tested two kinds of interpolations or
schemes

  1. 1.

A scheme that uses a QUICK third order interpolant in the bulk, away from the
interface and a simple first order upwind flux near the interface. We call
this scheme QUICK-UW.

  2. 2.

A scheme that uses a Superbee slope limiter [20]. for the flux in the bulk and
a more complex Superbee limiter tuned to a shifted interpolation point near
the interface. We call naturally call this scheme “Superbee”.

We describe each scheme in turn, starting with the QUICK-UW. For positive
advecting velocity and in the bulk we have

$\phi_{0}=\frac{3}{4}\phi_{-1/2}+\frac{3}{8}\phi_{1/2}-\frac{1}{8}\phi_{-3/2}$
(65)

while near the interface $\phi_{0}=\phi_{-1/2}$. For negative (to the left)
advecting velocity and in the bulk we have

$\phi_{0}=\frac{3}{4}\phi_{1/2}+\frac{3}{8}\phi_{-1/2}-\frac{1}{8}\phi_{3/2}$
(66)

while near the interface $\phi_{0}=\phi_{1/2}$ (see also [12])) In the
Superbee case, we start by defining a general family of interpolants that can
be expressed as

$f(\phi_{-3/2},\phi_{-1/2},\phi_{1/2},\phi_{3/2},{\rm
sign}(u_{f}))=\phi_{-1/2}+Sh/{2}$ (67)

for a positive advecting velocity where the slope $S$ is given by a slope-
limiter function $g$ such that $S=g(\phi_{-3/2},\phi_{-1/2},\phi_{1/2})$ while
for a negative advecting velocity

$f(\phi_{-3/2},\phi_{-1/2},\phi_{1/2},\phi_{3/2},{\rm
sign}(u_{f}))=\phi_{1/2}-Sh/{2}$ (68)

where $s=g(\phi_{-1/2},\phi_{1/2},\phi_{3/2})$.

To define the slope $S$ used in Superbee consider the factor $\alpha$ to be
the smallest in absolute value of

$\displaystyle\alpha^{+}=(\phi_{1/2,j}-\phi_{-1/2,j})/h,$
$\displaystyle\alpha^{-}=2(\phi_{-1/2,j}-\phi_{-3/2,j})/h.$ (69)

and the factor $\beta$ to be the smallest in absolute value of

$\displaystyle\beta^{+}=2(\phi_{1/2,j}-\phi_{-1/2,j})/h,$
$\displaystyle\beta^{-}=(\phi_{-1/2,j}-\phi_{-3/2,j})/h.$ (70)

Then one prescribes $S$ as the largest, in absolute value, of $\alpha$ and
$\beta$.

Figure 6: (a) The advected and advecting velocities in expression (58) are
both located on or near face $f$ (in this case $f=1-$). The reference control
volume $\Omega$ for the advected velocity component $u_{q}$ is also shown. The
arrangement is the same whatever the component $q$ but here the example of a
horizontal advection ( case $f=1-$ ) is taken. The value of the advected
velocity component $\bar{u}_{q}=\phi_{0}$ on the center of face $f$ (full
disk) is interpolated (see text) from the values $\phi_{p}$ of $u_{q}$ on the
nodes (open disks). (b) A more sophisticated interpolation predicts the value
$\phi(\hat{x})$ where $\hat{x}$ is at the center of the “donating” region
$\Omega_{D}$ (see text).

A slightly different estimate of the Superbee advected velocity is used near
the interface. First we extend the definition of the interpolants as we shall
predict $\bar{u}_{q}$ at a point $\hat{x}$ slightly upwind from $x_{f}$ using
a new function $\hat{f}$ so that

$\phi_{0}=\hat{f}(\hat{x},\phi_{-3/2},\phi_{-1/2},\phi_{1/2},\phi_{3/2},{\rm
sign}(u_{f})).$ (71)

We take $\hat{x}=x_{f}-u_{f}\tau/2$ to be the midpoint of the fluxed region.
(This is the region from which flow lines crossing the face originate see
Figure 6(b).) The extended interpolant is defined for positive velocity as

$\hat{f}(\hat{x},\phi_{-3/2},\phi_{-1/2},\phi_{1/2},\phi_{3/2},{\rm
sign}(u_{f}))=\phi_{-1/2}+S|\hat{x}-x_{-1/2}|$ (72)

and for negative velocity as

$\hat{f}(\hat{x},\phi_{-3/2},\phi_{-1/2},\phi_{1/2},\phi_{3/2},{\rm
sign}(u_{f}))=\phi_{1/2}-S|\hat{x}-x_{1/2}|$ (73)

The rationale behind this choice is as follows. For a time-independent
advecting velocity field $u({\bf x},t)$, the integrals in expression (47) can
be simplified

$\phi_{0}=\bar{u}_{q}=\frac{\int_{\Omega_{D}}u_{f}u_{q}\,{\rm d}{\bf
x}}{\int_{\Omega_{D}}u_{f}\,{\rm d}{\bf x}}.$ (74)

where ${\Omega_{D}}$ is the “donating region” (see [12]) from which flow lines
crossing the face originate. Then since we approximate the advecting velocity
$u_{f}$ by its midpoint value the integral can be further simplified as

$\bar{u}_{q}=\frac{1}{|\Omega_{D}|}\int_{\Omega_{D}}u_{q}\,{\rm d}{\bf x}$
(75)

Since the “midpoint” at $\hat{x}$ is the center of mass of the donating region
${\Omega_{D}}$ the interpolation expression (71) follows.

####  3.4.6 VOF-consistent Momentum advection: staggered grids.

In order to apply the above method on the staggered grid, the fluxes and
control volumes of the momentum must be in correspondence with the fluxes and
control volumes of the VOF color function. This is realized by recreating in
each velocity control volume the necessary color fraction data.

At the start of the velocity advection operations summarized by the operator
${\cal L}^{h}_{\rm conv}$ each velocity control volume say $\Omega_{i+1,j,k}$
overlaps two VOF control volumes $\Omega_{i,j,k}$ and $\Omega_{i+1,j,k}$ ( see
figure 7) The following operations are performed:

  1. 1.

Reconstruction of the volume fractions and density $C^{n}$ and $\rho^{n}$ in
the staggered control volumes $\Omega_{i+1/2,j,k}$ corresponding to the
momentum component $\rho u_{1}$, thus obtaining the “half fractions”
$\tilde{C}^{n}$ and $\tilde{\rho}^{n}$.

  2. 2.

Computation of the initial momentum component
$\tilde{g}_{1}^{n}=\tilde{\rho}^{n}u_{1}^{n}$.

  3. 3.

Split advection of the momentum in the $x$ direction using the method in the
previous sections to obtain the first component of the new momentum
$\hat{g}_{1}^{n+1}$.

  4. 4.

Split advection in the $x$ direction of the volume fractions and density in
the staggered control volumes using the volume of fluid method to obtain
predicted densities $\hat{\rho}^{n+1,1}$.

  5. 5.

Repeat the previous split advection operations for momentum and density in $y$
and $z$. At each time step, the sequence $x,y,z$ is permuted. Eventually one
obtains $\hat{g}_{1}^{n+1}$ and $\hat{\rho}^{n+1}$.

  6. 6.

Repeat the previous operations for the two other momentum components $g_{2}$
and $g_{3}$.

  7. 7.

Extraction of the new velocities components ${\bf u}^{*}=\hat{\bf
g}^{n+1}/\hat{\rho}^{n+1}$.

  8. 8.

In parallel, computation of $C^{n+1}$ from $C^{n}$ on the centered cells using
the VOF method.

We now cover the eight steps in detail. The reconstruction of the volume
fractions in the staggered cells (step 1) implies considering each control
volume for a velocity component, say $\Omega_{1}=\Omega_{i+1/2,j,k}$ for the
$u_{1}$ component, and finding the “half” volume fraction in each intersection
$\Omega_{i+1/2,j,k}\cap\Omega_{i,j,k}$ and
$\Omega_{i+1/2,j,k}\cap\Omega_{i+1,j,k}$ (see Figure 7). The “half-fractions”
are found using the reconstruction parameters $\alpha$, ${\bf m}$ in equation
(16) in the original centered control volume, and then computing the
corresponding volume fractions in the half cells. Addition of the two half
fractions leads to an estimate of the volume fraction
$\tilde{C}^{n}_{i+1/2,j,k}$ in the staggered cells, and thus of the densities
$\tilde{\rho}^{n}_{i+1/2,j,k}$.

In step 2, now that the staggered cell densities at the beginning of the time
step are known, obtain the momentum component
$g^{n}_{1,i+1/2,j,k}=\rho^{n}_{i+1/2,j,k}u^{n}_{1,i+1/2,j,k}$. This momentum
component can be advected (step 3) using the method described in the previous
section, to obtain $\hat{g}^{n+1}_{1,i+1/2,j,k}$ in a conservative manner
following equation (55). Similarly the volume fractions in the staggered cells
are also reconstructed and advected in direction $x_{1}$ in step 4. Then in
step 5 the other directionally split advections are performed. After each
direction, the newly obtained densities and advected velocities are used for
the next direction but the advecting velocities remain the same. The central
velocity stepping is discussed below equation (50) above. Finally the
advections are performed in each of the three systems of staggered cells (step
6) to obtain $\hat{C}^{n+1}_{i+1/2,j,k}=\hat{C}^{n,4}_{i+1/2,j,k}$ where the
notation $C^{n,m}$ was introduced in equation (20) and $m$ indicates which
direction of the split time step was performed. Then
$\hat{\rho}^{n+1}_{i+1/2,j,k}=\rho_{1}\hat{C}_{i+1/2,j,k}+\rho_{2}(1-\hat{C}_{i+1/2,j,k})$
obtains. Finally the ${\cal L}_{\rm conv}$ operator is obtained as (step 7)

${\cal L}_{\rm conv}=\frac{1}{\tau}(\hat{\bf g}^{n+1}-\tilde{\rho}^{n}{\bf
u}^{n})$ (76)

where as above we use the $\tilde{\rho}$ notation to indicate that the density
has been reconstructed in the staggered cells. In the absence of the other
operators in ${\cal L}_{2}$ this would result in

${\bf u}^{*}=\hat{\bf g}^{n+1}/\hat{\rho}^{n+1}.$ (77)

However (leaving aside the pressure and other contributions in equation (13))
the momentum on the staggered cells at the beginning of the next time step is
not equal to the momentum $\hat{\rho}^{n+1}{\bf u}^{n+1}$ at the end of the
previous time step but is rather defined as $\tilde{\rho}^{n+1}{\bf u}^{n+1}$
where $\tilde{\rho}^{n+1}$ is obtained at beginning of the time step from the
half-fractions. That is, $C^{n+1}$ is computed directly from $C^{n}$ on the
centered cells (step 8) and this allows to find $\tilde{C}^{n+1}$ and from it,
$\tilde{\rho}^{n+1}$ , by running again step 1. These densities
$\tilde{\rho}^{n+1}$ are not the same as the $\hat{\rho}^{n+1}$ densities of
the previous time step. This implies non conservation of the momentum. We note
that attempting to always use only the three sets $C^{n}_{i+1/2,j,k}$,
$C^{n}_{i,j+1/2,k}$, $C^{n}_{i,j,k+1/2}$ and evolve them by the VOF method on
the staggered cells would maintain conservation but result in the three
staggered grids evolving independently of each other and eventually diverging.

Figure 7: Computation of the staggered fractions from the half-fractions

###  3.5 Description of the other time-split terms

The other time-split terms in equation (12) and in the projection step (13)
are solved in a standard centered way. The density on the faces of the central
cells $\Omega_{i,j,k}$ is estimated using a centered average
$\rho_{i+1/2,j,k}=(\rho_{i+1/2,j,k}+\rho_{i+1/2,j,k})/2$. Although this is
less accurate and consistent than the usage of the half fractions
$\hat{\rho}_{i+1/2,j,k}$ described above the average is used in these terms
both for simplicity reasons and because tests have shown the usage of
$\hat{\rho}_{i+1/2,j,k}$ to lead to less stable simulations than the centered
average.

The velocities in the diffusion term are introduced in an explicit way.
Although this requires small time steps of the order $\rho h^{2}/\mu$ the
capillary restriction on time steps is usually even smaller, being of order
$\tau=(\rho h^{3}/\sigma)^{1/2}$. The two restrictions become of the same
order when $h\sim l_{\mu \sigma}$ where $l_{\mu\sigma}=\mu^{2}/(\sigma\rho)$
is the length at which the viscous and capillary terms balance. For water,
this length is of the order of 10 nanometers, and grids of that size are not
used in the flows we consider. However, should the velocities be treated in an
implicit manner, we do not believe this would change the conclusions of this
paper.

Surface tension is computed using the Continuous Surface Force method proposed
by [21], together with an estimate of the curvature through the computation of
height functions, in a manner that closely follows the method of [22]. The
external forces in equation (12) are only gravity and are computed in a
trivial manner with $\frac{1}{\rho^{n+1}}{\cal L}_{ext}={\bf g}$, where
gravity ${\bf g}$ is a constant.

##  4 Testing and Validation

###  4.1 Consistent cylinder advection

An elementary test of our method, that mostly verifies that the coding has
been performed correctly, considers a uniform planar velocity field
$u_{1}=u_{2}=1.6\times 10^{-2}$ and a droplet of density $\rho_{l}=10^{9}$ in
gas at density $\rho_{g}=1$ with a CFL number of $0.0256\sqrt{2}$. Viscosity
and surface tension are set to zero in this test. The diameter of the droplet
is $D/h=3.2$. The unit domain is discretized on a $16\times 16$ grid. The
droplet shapes that result are shown on the left of Figure 8. The
irregularities seen in the advected droplet are due to the roughness of the
VOF approximation at such low resolutions. We repeat the test with conditions
close to air/water: viscosity $\mu_{l}=0.1$, $\mu_{g}=0.002$ and $\rho_{g}=1$,
$\rho_{l}=10^{3}$ with identical results: a viscosity contrast will not
generate numerical instabilities on a uniform velocity field, as shown on the
right of Figure 8.

Figure 8: Large-density-ratio droplet in a uniform velocity field: a droplet
with $D/h=3.2$ grid points per diameter is advected in 3D in the plane $z=0$
(see text); left: density ratio $10^{9}$ without viscosity, right: density
ratio $10^{3}$ with viscosity

###  4.2 Sudden acceleration of a cylinder at large density contrast

A test that is often included in studies of momentum-conserving methods [2, 3,
4, 5, 6] and other methods designed to improve the stability of two-phase flow
computations [23] is to initialize a droplet of very high density at velocity
${\bf u}_{l}({\bf x})=U_{0}{\bf e}_{x}$ with the other lighter fluid at rest,
so that ${\bf u}_{g}({\bf x})=0$. Surface tension and viscosity are not
present as in the previous test, the only difference being the discontinuity
of the initial velocity on the interface. The initial velocity condition
amounts to a vortex sheet on the surface of the cylinder. After the first time
step, the projection method (13) adds a dipole potential flow so that ${\bf
u}_{g}=U_{0}{\bf e}_{x}+\tau\nabla p/\rho_{g}$ in the gas around the droplet,
identical to the dipole flow around a solid object. However in addition to the
dipole flow there could be small multipole components to the outer gas flow.
There are three ways to estimate the perturbations created on the interface.

Figure 9: Deformation of a droplet moving through a light fluid. From left to
right and top to bottom: $r=10,10^{3},10^{6},10^{9}$. Black shape at $t=0$,
blue at $t=D/(2U)$ and red at $t=D/U$. The grid resolution is $D/h=20$. Figure
10: Deformation of a droplet moving through a light fluid at $r=10^{9}$. The
grid resolution is $D/h=20$.

A vortex sheet is unstable with respect to the Kelvin-Helmholtz instability.
The main results about the amplitude of the instability are as follows. Let
$A(t)$ be the maximum deviation of the interface which has a dimension of
length. For flows with no viscosity and surface tension as is the case here,
the Kelvin-Helmholtz instability amplitude $A(t)$ should grow exponentially
from a perturbation of wavenumber $k$ as $A(t)\sim A(0)\exp(s_{k}t)$ with
$s_{k}=|\Delta U|k/\sqrt{r}$ where $\Delta U$ is the tangential velocity
difference and $r=\rho_{l}/\rho_{g}$. In the limit of large $r$ the growth
rate becomes small. Since the maximum wavenumber on the grid is $\pi/h$ an
estimate of the growth rate of the small wavelength instabilities is $\pi
U_{0}N/(D\sqrt{r})$ where $N$ is the number of grid points per diameter. After
advection by a droplet diameter, the elapsed time is $\Delta t=D/U_{0}$. For
typical values in the literature of $r=10^{6}$ and an arbitrary value of
$N=32$ the amplitude growth would be

$\exp(s_{\rm max}\Delta t)=\exp\left(\frac{\pi
N}{\sqrt{r}}\right)=\exp(0.032\pi)=1.1058\ldots$ (78)

which means the amplitude should grow by 10% after advection by one diameter
and by $e^{1/2}$ after the typical advection by 5 diameters.

Beyond the linear growth stage of the Kelvin-Helmholtz instability, there is a
self-similar, non-linear growth stage for which dimensional analysis implies
that $A(t)\sim|\Delta U|t/\sqrt{r}$ [24]. By this argument also the
perturbation of the cylinder should remain of order $A(\Delta t)\sim
D/\sqrt{r}$ after advection by a droplet diameter. One should note that the
self-similar growth is obtained for vanishing boundary layer thickness but we
precisely expect such a velocity discontinuity in the present case with no
viscosity.

A final, physical deformation is expected from the spatial pressure variation
induced by the dipole gas flow. This variation involves a larger pressure at
the aft and fore stagnation points and a lower pressure along the “equator” of
the droplet [25]. The resulting integrated stress is of order $\rho_{g}U^{2}$
resulting in the growth of the droplet deformation as $A(t)\sim
U^{2}t^{2}/(Dr)$ and after advection by a droplet diameter as $A(\Delta t)\sim
D/r$. This growth is observed experimentally [26] and results in an
elliptically shaped drop, albeit of much smaller amplitude than the two former
Kelvin-Helmholtz-related growth mechanisms.

The results are shown on Figure 9 for several times and density ratios in a
manner comparable to [2]. For these numerical experiments, the exact manner of
initializing the velocity fields has some importance. It is important not to
allow in the first time step some gas velocity in the liquid to avoid large
pressure gradients in the liquid. Thus the density is initialized to
$\rho_{l}$ to machine accuracy using the Vofi library [13, 14] within a disk
implicitly defined by $x^{2}+y^{2}+z^{2}<R^{2}$ and the velocity is
initialized to 1 for all the velocity nodes inside a disk implicitly defined
by $x^{2}+y^{2}+z^{2}<(R+nh)^{2}$, where $n$ is the size of the “halo” in
number of grid points. The velocity in the other nodes is initialized to 0.
The tests shown were performed with $n=1$. Increasing the size of the halo
from $n=0$ to $n=1$ improves the results in the first phase of the droplet
motion, but not in the late phase.

The velocity of the motion has been oriented on the diagonal as in [2]. The WY
scheme is used with a QUICK interpolant. The droplet deforms little after
advection by one droplet diameter (Figure 9). For longer advection the
deformation is worse but, as explained above, this is to be expected at high
resolution, except in the $r=10^{9}$ case for which we show advection by 5
diameters in Figure 10. If the non-momentum-conserving method is used, the
droplet deforms rapidly and the simulation breaks down. As expected the
results are better than those of [23], that uses the skew-symmetric scheme but
without momentum-conserving-like schemes, but worse than those of [2, 3, 4,
5]. It is also seen, as expected from the developments above, but perhaps
contrary to intuition, that it is easier to have larger density contrasts.
Increasing the number of grid points also helps.

A comment can be made on the nature of the instabilities seen. If the gas
velocity numerically diffuses inside the liquid, then some vorticity may
penetrate into the liquid despite the fact that in inviscid flow vorticity
should remain confined on the interface. If this happens, the growth rate of a
single-phase Kelvin-Helmholtz instability inside the liquid is the much larger
value $\widetilde{s}_{k}=|\Delta U|k/2$, without the $\sqrt{r}$ factor at
denominator. (The derivation of this fact may be found for example in [27].)
The instability seen for the $r=10^{9}$ case for the long advection case in
Figure 10 is a mode 3, which excludes the large linear growth discussed above
for small wavelenghts rather points to the velocity diffusion mechanism.
Minimizing numerical diffusion should then be an important quality for this
test, while on the other hand numerical diffusion may stabilize other
difficult test cases/situations.

###  4.3 Sheared layer

In order to better analyze the behavior of the methods in flow under shear,
such as vortex sheets and Kelvin Helmholtz instabilities, we setup directly a
planar, parallel shear flow in a $(-1/2,1/2)^{2}$ domain with the following
initial condition

$\displaystyle u_{1}$ $\displaystyle=15$ $\displaystyle\qquad{\rm
if\;\;}|x_{2}|>1/10$ $\displaystyle u_{1}$ $\displaystyle=1$
$\displaystyle\qquad{\rm if\;\;}|x_{2}|<1/10$ $u_{2}=0.01\,\sin(2\pi
x_{1})\exp(-20x_{2}^{2})$

and $\rho=1$ if $|x_{2}|>1/10$ and $\rho=10^{3}$ otherwise. This flow is
similar to a liquid sheet in high velocity gas. The flow is simulated until
time $t_{f}=2$ unless the simulation blows up at an earlier time
$t_{b}<t_{f}$. Table 1 shows the numerical schemes that have been used in a
few simulations together with the fraction $t_{b}/t_{f}$ of the final time
that has been reached. All simulations have been performed on a $128\times
128$ grid with a CFL of 0.03 and a density ratio of 1000. The new method (MC)
is systematically more stable than the standard method for the two
combinations shown in Table 1: a) the WY VOF scheme and the QUICK-UW velocity
interpolation schemes and b) the CIAM VOF with the Superbee slope limiters.
The state of the two simulations with the first combination of schemes is
shown in Figure 11, just before breakdown at time $t_{b}=0.34$ (or
$t_{b}/t_{f}=0.17$) for the (MC) method, and at time $t_{b}=0.16$ (or
$t_{b}/t_{f}=0.08$) for the standard method. We have also tested a number of
other combinations, for example WY & Superbee that turns out to be very
unstable. We have also performed simulations with smaller $64\times 64$ and
$32\times 32$ grids that yield similar results. Finally, we note that this
case is also unstable in a single-phase configuration when using the QUICK
third-order velocity interpolation.

Method | VOF Scheme | velocity scheme | fraction completed
---|---|---|---
standard | WY | QUICK-UW | 0.08
MC | WY | QUICK-UW | 0.17
standard | CIAM | Superbee | 0.50
MC | CIAM | Superbee | 1
Table 1: Percentage of completion of the shear layer test with various
methods. All simulations are performed on a $128\times 128$ grid with a CFL of
0.03 and a density ratio of 1000. Figure 11: The state of the simulation of
the sheared layer just before eventual breakdown of the simulation with the
combination of WY VOF scheme and QUICK-UW interpolation. Left: (MC) method,
right: standard method.

###  4.4 Falling Rain Drop

Similar to other tests of large-density-contrast schemes, we study the fall of
a three-millimeter water drop in air (in other words, a relatively large
raindrop). The terminal velocity of such a raindrop is about $U=8m/s$ (see
ref. [28]) and thus the Weber number is about ${\rm We}\,=\rho_{\rm
air}U^{2}d/\sigma\simeq 3.2$ and the Reynolds number of the air flow is about
${\rm Re}\,=\rho_{\rm air}Ud/\mu_{\rm air}\simeq 1455$. The problem setup is
shown on Fig. 12 and the parameters are given in Table 2. A time-dependent
inflow velocity is set, with the velocity $U_{0}(t)$ being adjusted by a small
controller program that eventually stabilizes the droplet at mid-height of the
domain.

$\rho_{\rm air}$ | $\rho_{\rm water}$ | $\mu_{\rm air}$ | $\mu_{\rm water}$ | $\sigma$ | $d$ | $g$
---|---|---|---|---|---|---
$\left(kg/m^{3}\right)$ | $\left(kg/m^{3}\right)$ | $\left(Pa\,s\right)$ | $\left(Pa\,s\right)$ | $\left(N/m\right)$ | $(m)$ | $(m/s^{2})$
1.2 | $0.9982\,10^{3}$ | $1.98\,10^{-5}$ | $8.9\,10^{-4}$ | $0.0728$ | $3\,10^{-3}$ | $9.81$
Table 2: Parameter values used in the simulation of a falling water droplet in
air.

For such a Weber number the capillary forces dominate and the droplet should
remain approximately spherical. Instead at moderate resolution ($D/h=15$ to 60
grid points per diameter were tested) and without the momentum-conserving
scheme described in this paper the droplet deforms catastrophically as shown
on Fig. 13.

Figure 12: Problem setup for falling rain drop test case. Boundary conditions
at the top and bottom are a uniform inflow and outflow velocity $U_{0}(t)$.
Boundary conditions on the side are free slip (no shear stress). Figure 13:
Rain drop test case: catastrophic breakup with non-conserving formulation,
$D/h=30$.

This effect was already found by [11] in a similar case, the sudden
interaction of a droplet at rest with uniform gas flow. For a Weber number
also around 3, the droplet should remain near spherical but ref. [11] reports
a similar catastrophic deformation and gives the following explanation. To
start with, we neglect gravity and viscous effects at this relatively large
Reynolds number. Also, we are interested in steady-state flow. On the axis and
near the hyperbolic stagnation point at the front of the droplet one has
$u_{2}=0$ for the transverse (radial) velocity and for the axial momentum
balance

$u_{1}\partial_{1}u_{1}=-\frac{1}{\rho}\partial_{1}p.$ (79)

Because of the low viscosity and large density ratio, it is not possible for
the air flow to immediately entrain the water, so the fluid velocity is
significantly smaller in the water. In the air the acceleration near the
stagnation point is $U^{2}/D$ and the pressure gradient is

$\partial_{1}p\sim\rho_{a}U^{2}/D.$ (80)

The pressure gradient in the water is much smaller, however if in a mixed cell
the water density multiplies the air acceleration $U^{2}/D$, so that

$\partial_{1}p\sim\rho_{w}U^{2}/D,$ (81)

then a large pressure gradient results in the mixed cell or cells. This large
pressure gradient results in a large pressure inside the droplet near the
front stagnation point, as shown in Figure 14. This large pressure is balanced
by surface tension only for a sufficiently large curvature near the droplet
front. This explains the presence of a “dimple” often seen in low resolution
simulations of the falling drop.

|
---|---
(a) | (b)
Figure 14: The origin of the pressure peak in the front of the droplet. (a) The profile of the pressure on the axis a few timesteps after initialisation with the standard, non-momentum conserving method (red) and the present method (blue). (b) The pressure distribution immediately after the start of the simulation using the standard, non-momentum-conserving method. The pressure peak does not result yet in the formation of a dimple. In all figures $D/h=16$. |
---|---
(a) | (b)
Figure 15:  (a) The pressure profile on the axis a few timesteps after
initialisation with the present method at various resolutions: $D/h=8$ (red) ,
$16$ (blue) and $32$ (green). (b) The pressure distribution immediately after
the start of the simulation using the present method and $D/h=16$.

Applying the method described in this paper brings a considerable improvement,
as shown by a comparison of Figures 14 and 15. Simulations fall in three
categories: those that blow up anyway, those that have a marked peak in
kinetic energy as a function of time, associated with deformed interface
shapes, and those that keep physical values of the kinetic energy and smooth
shapes. The most stable combinations of schemes are CIAM advection with
superbee limiter and WY advection with QUICK limiter. The CIAM advection with
the superbee limiter appears also to be very diffusive.

|
---|---
(a) | (b)
Figure 16: Flow field around the 3 mm droplet with 60 grid points per
diameter. (a) The velocity magnitude. It is seen that even at this highest
resolution there are only three points in the boundary layer. (b) The
vorticity magnitude. The marked separation of the boundary layers is observed
with a more complex vortical region in the wake.

Visualisation of the flow around the droplets (Figure 16) contributes to
explain why it is so challenging. It can be seen that the boundary layers are
very thin, questioning our approximation that fluid velocity is continuous
across the interface.

In order to validate our simulation, we study the convergence upon grid
refinement. Three grids are used, with $D/h=15,\,30$ and $60$. A $12$-mm cubic
box is used with the $3$-mm droplet at the center. We study the convergence of
the terminal velocity and that of the shape. For the latter, we use as a
descriptor of the shape the three moments of inertia $I_{m}$ defined by

$I_{m}=\int_{\cal D}Hx_{m}^{2}{\rm d}{\bf x}\;,\quad 1\leq m\leq 3,$ (82)

where ${\cal D}$ is the domain used for the computation and $x_{m}$ is
relative to the center of mass. The convergence of the moments of inertia and
terminal velocity is shown on Figure 17. The velocity seems to converge to a
value around $7\,m/s$, to be compared with a value of $8.06\,m/s$ found by the
authors of ref. [28]. This is not surprising since the experimental terminal
velocity would be modified for a droplet constrained in a finite-size box.

|
---|---
(a) | (b)
Figure 17: Convergence of simulations. (a) Evolution of the terminal velocity
with grid refinement. (b) Evolution of the three moments of inertia with grid
refinement. Figure 18: Flattening of the droplet with increasing equivalent
diameter (see text). From left to right $D_{e}=3,\,4.6,\,6.4$ and $8\,mm$.

From the values of the moments of inertia, the horizontal and vertical extents
of the droplet (respectively $D_{r}$ and $D_{x}$) can be found and compared to
the values found by the authors of [29]. We find $D_{r}=3.1\,mm$ and
$D_{x}=2.6\,mm$, while ref. [29] concludes that drops are quasi-spherical for
an equivalent diameter $D_{e}\leq l_{c}$ and deformed for $D_{e}>l_{c}$, where
$l_{c}=(\sigma/\rho g)^{1/2}$ is the capillary length, $l_{c}=2.7\,mm$ for
water. Indeed repeating the simulations for larger drops we find increased
flattening as shown on Figure 18.

###  4.5 Atomizing air and water planar jets

We also test the capability of the VOF-consistent momentum-conserving scheme
to simulate complex air-water flows with an unstable shear flow. For that
purpose we repeat the setup of reference [30]. Two jets of air and water are
entering the computational domain from the left of Figure 19 at velocities
comparable to those of experiments. However in order to save computational
time the domain is smaller than in experiments. Physical properties of air and
water are identical to those of the falling raindrop case given in Table 2.
The flow and domain characteristics are given in Table 3 including a gas
boundary layer and separator-plate size identical to those of reference [30].
The notations are as in reference [30]: $H_{p}$ is the thickness of the jet of
phase $p$, there is a separator plate of thickness $l_{y}$ and a gas boundary
layer of thickness $\delta_{g}$. The two streams have equal thickness
$H_{l}=H_{g}$. The dimensionless parameters are given in Table 4. The CIAM
advection method has been used. The number of grid points in the layer
$H_{l}/h=16$ is relatively small (compared to the $H_{l}/h=32$ in the smallest
(coarsest) simulation of reference [30]). It is thus all the more remarkable
that the simulation is stable since using a smaller number of grid points
usually increases the trend towards instability. It is interesting to note
that the VOF calculations accounts for 31.5% of the total time, while the
inversion of the Poisson operator for the pressure accounted for 51.5%. The
whole simulations runs overnight on a present-day workstation.

Figure 19: Atomizing layer with air/water properties. $U_{l}$ | $U_{g}$ | $H_{l}$ | $h$ | $l_{y}$ | ${\delta_{g}}/{l_{y}}$
---|---|---|---|---|---
$(m/s)$ | $(m/s)$ | $(m)$ | $(m)$ | $(m)$ | $(-)$
$1$ | $25$ | $4\,10^{-3}$ | $2.5\,10^{-4}$ | $2.5\,10^{-4}$ | $2$
Table 3: Physical parameters (defined in the text) for the atomizing layer setup. The fluid properties are the same as in Table 2.  $M$ | $r$ | $m$ | ${\rm Re}\,_{g,\delta}$ | ${\rm We}\,_{g,\delta}$ | ${\rm Re}\,_{g}$
---|---|---|---|---|---
$\rho_{g}U_{g}^{2}/\left(\rho_{l}U_{l}^{2}\right)$ | $\rho_{l}/\rho_{g}$ | $\mu_{l}/\mu_{g}$ | $\rho_{g}U_{g}H_{g}/\mu_{g}$ | $\rho_{l}U_{l}H_{l}/\mu_{l}$ | $\rho_{g}U_{g}^{2}H_{g}/\sigma$
$0.75$ | $831.8$ | $45$ | $757.6$ | $5.151$ | $6061$
Table 4: Dimensionless parameters for the atomizing layer setup.

##  5 Conclusion

We have presented and tested a simulation method for multiphase flow that
shows increased stability at large density contrasts and large Reynolds
numbers. The method is closely related to some variants of the VOF method
including an implementation of the WY advection method. The increased
stability is evidenced in the test case of a single $3$-mm droplet of water
falling in air, a typical raindrop. It is a reflection on the challenging
nature of multiphase flow that such complex methods apparently need to be
implemented to resolve such an everyday and simple phenomenon.

The method comes with a significant saving of computer time, since for similar
problems with raindrops, our attempts with a non-momentum-conserving VOF
approach led to catastrophic deformation of the drop or strong dimple
formation. These problems have also been observed by us using other non-VOF-
consistent and non-momentum-conserving methods such as the one of [22]. In
that case whenever less than 200 grid points per diameter are used numerically
stable air-water drops accelerated at moderate Weber number cannot be found.
However for higher resolutions they can be computed without difficulty as also
found by the authors of ref. [31]. Here, approximate solutions accurate within
15% are found with only 15 points per diameter, and non-divergent computations
are found with as little as 2 points per diameter.

A particular advantage of the method is that it is conserving mass at the
accuracy at which discrete incompressibility is enforced and opens a
perspective for similar momentum conservation using WY advection. The method
nevertheless is more complex and costly than a colocated method. This opens
the perspective for systematic development of other methods with different
grid arrangements. Another perspective is the potential of stable methods with
large density contrasts, exact mass and momentum conservation and small
droplets, that could be smoothly merged into models that represent the small
droplets as Lagrangian Point Particles [32].

##  6 Acknowledgements

This work has been supported by the ANR MODEMI project (ANR-11-MONU-0011)
program and grant SU-17-R-PER-26-MULTIBRANCH from Sorbonne Université.

This work was granted access to the HPC resources of TGCC-CURIE, TGCC-IRENE
and CINES-Occigen under the allocations t20152b7325, t20162b7760, 2017tgcc0080
and A0032B07760, made by GENCI and TGCC. The authors would also acknowledge
the MESU computing facilities of Sorbonne Université.

We would thank Dr. W. Aniszewski, Dr S. Dabiri, Dr Jiacai Lu and Dr. P. Yecko
for their contribution to the development of the code _PARIS-Simulator_ , and
we thank Dr. W. Aniszewski, Dr. V. Le Chenadec, Mr. C. Pairetti, Dr. S.
Popinet and Dr. S. Vincent for useful conversations on the topics of this
paper.

Finally, the simulation data are visualized by the software VisIt developed by
the Lawrence Livermore National Laboratory.

## References

## References

  * [1] M. Rudman, A volume-tracking method for incompressible multi-fluid flows with large density variations, Int. J. Numer. Meth. Fluids 28 (1998) 357–378.
  * [2] M. Bussmann, D. B. Kothe, J. M. Sicilian, Modeling high density ratio incompressible interfacial flows, in: ASME 2002 Joint US-European Fluids Engineering Division Conference, American Society of Mechanical Engineers, 2002, pp. 707–713.
  * [3] O. Desjardins, V. Moureau, Methods for multiphase flows with high density ratio, Center for Turbulent Research, Summer Programm 2010 (2010) 313–322.
  * [4] M. Raessi, H. Pitsch, Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method, Computers & Fluids 63 (2012) 70–81.
  * [5] V. Le Chenadec, H. Pitsch, A monotonicity preserving conservative sharp interface flow solver for high density ratio two-phase flows, J. Comput. Phys. 249 (2013) 185–203.
  * [6] G. Vaudor, T. Menard, W. Aniszewski, M. Doring, A. Berlemont, A consistent mass and momentum flux computation method for two phase flows. Application to atomization process, Computers & Fluids 152 (2017) 204–216.
  * [7] R. Saurel, R. Abgrall, A multiphase Godunov method for compressible multifluid and multiphase flows, J. Comput. Phys. 150 (2) (1999) 425–467.
  * [8] G. Allaire, S. Clerc, S. Kokh, A five-equation model for the simulation of interfaces between compressible fluids, J. Comput. Phys. 181 (2) (2002) 577–616.
  * [9] M. Pelanti, K. M. Shyue, A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves, J. Comput. Phys. 259 (2014) 331–357.
  * [10] F. Xiao, Large eddy simulation of liquid jet primary breakup, Ph.D. thesis, Loughborough University (2012).
  * [11] F. Xiao, M. Dianat, J. J. McGuirk, Large eddy simulation of single droplet and liquid jet primary breakup using a coupled level set/volume of fluid method, Atomisation and Sprays 24 (2014) 281–302.
  * [12] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct Numerical Simulations of Gas-Liquid Multiphase Flows, Cambridge University Press, 2011.
  * [13] S. Bnà, S. Manservisi, R. Scardovelli, P. Yecko, S. Zaleski, Numerical integration of implicit functions for the initialization of the VOF function, Computers & Fluids 113 (2015) 42–52.
  * [14] S. Bnà, S. Manservisi, R. Scardovelli, P. Yecko, S. Zaleski, Vofi — a library to initialize the volume fraction scalar field, Computer Physics Communications 200 (2016) 291–299.
  * [15] R. Scardovelli, S. Zaleski, Analytical relations connecting linear interfaces and volume fractions in rectangular grids, J. Comput. Phys. 164 (2000) 228–237.
  * [16] J. Li, Calcul d’interface affine par morceaux (piecewise linear interface calculation), C. R. Acad. Sci. Paris, série IIb, (Paris) 320 (1995) 391–396.
  * [17] R. Scardovelli, S. Zaleski, Interface reconstruction with least-square fit and split lagrangian-eulerian advection, Int J. Numer. Meth. Fluids 41 (2003) 251–274.
  * [18] G. D. Weymouth, D. K. P. Yue, Conservative Volume-of-Fluid method for free-surface simulations on Cartesian-grids, J. Comput. Phys. 229 (8) (2010) 2853–2865.
  * [19] S. Popinet, Basilisk, a Free-Software program for the solution of partial differential equations on adaptive Cartesian meshes (2018).
URL http://basilisk.fr

  * [20] P. L. Roe, Some contributions to the modelling of discontinuous flows, in: Large-scale computations in fluid mechanics, 1985, pp. 163–193.
  * [21] J. Brackbill, D. B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335–354.
  * [22] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, J. Comput. Phys. 228 (2009) 5838–5866.
  * [23] D. Fuster, An energy preserving formulation for the simulation of multiphase turbulent flows, J. Comput. Phys. 235 (2013) 114–128.
  * [24] R. Blumenthal, J. Hoepffner, S. Zaleski, Self-similar wave produced by local perturbation of the Kelvin-Helmholtz shear-layer instability, Phys. Rev. Lett 106 (2011) 104502.
  * [25] R. Clift, J. R. Grace, M. E. Weber, Bubbles, Drops and Particles, Academic Press, New York, 1978.
  * [26] L. Opfer, I. Roisman, J. Venzmer, M. Klostermann, C. Tropea, Droplet-air collision dynamics: Evolution of the film thickness, Phys. Rev. E 89 (1) (2014) 013023.
  * [27] J.-P. Matas, S. Marty, A. Cartellier, Experimental and analytical study of the shear instability of a gas-liquid mixing layer, Physics of fluids 23 (9) (2011) 094112.
  * [28] R. Gunn, G. D. Kinzer, The terminal velocity of fall for water droplets in stagnant air, Journal of Meteorology 6 (4) (1949) 243–248.
  * [29] É. Reyssat, F. Chevy, A. L. Biance, L. Petitjean, D. Quéré, Shape and instability of free-falling liquid globules, Europhys. Lett. 80 (3) (2007) 34005–6.
  * [30] Y. Ling, D. Fuster, S. Zaleski, G. Tryggvason, Spray formation in a quasiplanar gas–liquid mixing layer at moderate density ratios: A numerical closeup, Phys. Rev. Fluids 2 (1) (2017) 014005.
  * [31] M. Jain, R. S. Prakash, G. Tomar, R. V. Ravikrishna, Secondary breakup of a drop at moderate Weber numbers, in: Proc. R. Soc. A, Vol. 471, issue 2177, The Royal Society, 2015, p. 20140930.
  * [32] Y. Ling, S. Zaleski, R. Scardovelli, Multiscale simulation of atomization with small droplets represented by a lagrangian point-particle model, Int. J. Multiphase Flow 76 (2015) 122 – 143.  doi:http://dx.doi.org/10.1016/j.ijmultiphaseflow.2015.07.002.
URL http://www.sciencedirect.com/science/article/pii/S0301932215001524