# PArallel, Robust, Interface Simulator (PARIS)

W. Aniszewski  T. Arrufat  M. Crialesi-Esposito  S. Dabiri  D. Fuster  Y. Ling
J. Lu  L. Malan  S. Pal  R. Scardovelli  G. Tryggvason  P. Yecko  S. Zaleski
Sorbonne Université & CNRS, UMR 7190,
Institut Jean Le Rond d’Alembert, F-75005, Paris, France  CMT-Motores
Térmicos, Universitat Politécnica de Valéncia, Camino de Vera, s/n, Edificio
6D, Valencia, Spain  Mechanical Engineering, Purdue University, West
Lafayette, IN, USA  Mechanical Engineering, Baylor University, Waco, TX 76706,
USA  Mechanical Engineering, University of Cape Town, South Africa  Mechanical
Engineering, Johns Hopkins University, Baltimore, MD, USA  DIN - Lab. di
Montecuccolino, Università di Bologna, I-40136 Bologna, Italy  Cooper Union,
New York City, USA

###### Abstract

Paris (PArallel, Robust, Interface Simulator) is a finite volume code for
simulations of immiscible multifluid or multiphase flows. It is based on the
“one-fluid” formulation of the Navier-Stokes equations where different fluids
are treated as one material with variable properties, and surface tension is
added as a singular interface force. The fluid equations are solved on a
regular structured staggered grid using an explicit projection method with a
first-order or second-order time integration scheme. The interface separating
the different fluids is tracked by a Front-Tracking (FT) method, where the
interface is represented by connected marker points, or by a Volume-of-Fluid
(VOF) method, where the marker function is advected directly on the fixed
grid. Paris is written in Fortran95/2002 and parallelized using MPI and domain
decomposition. It is based on several earlier FT or VOF codes such as Ftc3D,
Surfer or Gerris. These codes and similar ones, as well as Paris, have been
used to simulate a wide range of multifluid and multiphase flows.

###### keywords:

Multiphase flows , Multifluid flows , Free-surface flows , Navier-Stokes
equations , Front Tracking , Volume of Fluid , Surface Tension

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: PArallel Robust Interface Simulator — Paris
Licensing provisions: GPLv3.
Programming language: Fortran95/2002. Parallelized using MPI and domain
decomposition.
Nature of problem:
Paris is a free code, or software, for computational fluid dynamics (CFD) of
multiphase flows (or computational multiphase fluid dynamics (CMFD)),
specialized in the numerical simulation of interfacial fluid flows, involving
droplets, bubbles and waves, as described in the book by Tryggvason,
Scardovelli and Zaleski [1]. It solves the Euler or Navier-Stokes equations in
the one-fluid formulation of two-phase flow, including a surface tension
force. It computes complex flows such as fast atomizing jets or droplets,
expanding cavitation bubble clusters and multiphase flow through porous media.
Solution method:
The code mostly implements the methods described in the book by Tryggvason,
Scardovelli and Zaleski [1]. Time stepping is performed using a first-order or
a second-order in time predictor-corrector method with an explicit projection
step for the pressure. Spatial discretization is by finite volumes on a
regular cuboid grid. Interface tracking is performed with the Front-Tracking
(FT) method or the Volume-of-Fluid (VOF) method. In the VOF version Paris uses
either the Lagrangian-Explicit (LE) advection method or the exactly mass-
conserving method of Weymouth and Yue [2]. The normal computation is performed
using the Mixed-Youngs-Centered (MYC) scheme. A mass-momentum advection method
has been also implemented that is consistent with the VOF advection [3].
Curvature is computed with the Height Function (HF) method. This is combined
with the balanced Continuous Surface Force (CSF) method to compute surface
tension forces. If the dynamics of a phase can be neglected, Paris can also
run as a free-surface code by specifying a homogeneous pressure, at most
varying with time, in the neglected phase. In the case of atomizing jets, an
algorithm has been implemented in Paris that can detect isolated droplets,
advects them as Lagrangian point-particles and possibly merge them again with
the main stream.
Additional comments:
Paris is extended from or inspired by the following codes:

  * 1.

Ftc3D: Front Tracking code for 3D simulations by Gretar Tryggvason and Sadegh
Dabiri

  * 2.

Surfer: VOF code for 3D simulations by Stephane Zaleski, Jie Li, Ruben
Scardovelli and others

  * 3.

Gerris: multiphase flow solver with Adaptive Mesh Refinement (AMR) by Stephane
Popinet

## References

  * [1] G. Tryggvason, R. Scardovelli, and S. Zaleski.  Direct Numerical Simulations of Gas-Liquid Multiphase Flows.  Cambridge University Press, 2011.
  * [2] G. D. Weymouth and Dick K. P. Yue.  Conservative Volume-of-Fluid method for free-surface simulations on Cartesian-grids.  Journal of Computational Physics, 229(8):2853–2865, April 2010\.
  * [3] T. Arrufat, M. Crialesi-Esposito, D. Fuster, Y. Ling, L. Malan, S. Pal, R. Scardovelli, G. Tryggvason, S. Zaleski, A mass-momentum consistent, Volume-of-Fluid method for incompressible flow on staggered grids,  Computers & Fluids, 215, 104785, 2021.

##  1 Introduction

Computations of the unsteady motion of multifluid flows, where two or more
immiscible fluids or thermodynamic phases flow while separated by sharp
interfaces, date back to the earliest days of computational fluid dynamics
(see [1, 2] for reviews). However, early simulations were restricted to
relatively small and idealized problems. As computer power has continued to
grow, it has been increasingly possible to conduct Direct Numerical
Simulations (DNS), defined as fully resolved and verified simulations of a
validated system of equations that include non-trivial length and time scales.

A few authors of Paris [3] have been involved in the development of free codes
for DNS of two-phase flows, such as Surfer [4], Gerris [5], dating back a
couple of decades, and more recent ones such as Basilisk [6]. The new project
Paris is a joint effort with the aim to illustrate most of the methods
described in [2] (the latter book will be denoted “TSZ” in what follows to
simplify citations and we will often refer to “TSZ” for developments about the
numerical methods) and it relies heavily on the previous codes Ftc3D and
Surfer. The software can run independently either with the Front-Tracking
method or the VOF method. The mass-momentum consistent advection method, the
free-surface option and the Lagrangian point-particles (LPP) algorithm have
been specifically designed for Paris and were not present in the previous
codes.

Furthermore, a software platform such as Paris can make it possible to explore
in the near future the development of combined models based on both Front-
Tracking and Volume-of-Fluid methods.

While DNS of multiphase flows are becoming increasingly common, most research
groups need to devote a considerable amount of time to code development. For
new groups, the need to develop a suitable simulation tool can be a
significant barrier to entry. Paris is a free code that is intended to be
relatively simple to use and modify even by beginner users, yet it has
sufficient capabilities to allow state-of-the-art studies of typical
multiphase flow problems.

##  2 Navier–Stokes equations with interfaces

###  2.1 Basic equations

In the three-dimensional space, the locus of the interface is a smooth surface
$S$ which separates the two fluid phases. Accordingly, we assume that the
interface is an object of zero thickness. This latter assumption constitutes
the “sharp interface” approximation. In this approximation, the phases are
implicitly located by a Heaviside function $\chi({\bf x},t)$ defined such that
fluid 1 corresponds to $\chi=1$ and fluid 2 to $\chi=0$. Viscosity and density
$\mu$ and $\rho$ are space and time dependent and given by

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

In the case with no phase change, mass conservation implies that the interface
advances at the speed of the flow, that is

$V_{S}={\bf u}({\bf x},t)\cdot{\bf n}$ (2)

where ${\bf u}({\bf x},t)$ is the local fluid velocity and ${\bf n}$ a unit
normal vector perpendicular to the interface. Equivalently, this condition on
the interface motion can be expressed, in weak form, as

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

which expresses the fact that the singularity of $\chi$, located on $S$, moves
at velocity $V_{S}={\bf u}\cdot{\bf n}$. We refer the reader to the
literature, in particular TSZ, for additional developments on interface
geometry.

For incompressible flows, which we will consider in what follows, we have

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

The Navier–Stokes equations for incompressible, Newtonian flow with surface
tension may conveniently be written in a conservative form, expressing the
momentum balance, or in a non conservative, Lagrangian form. The first form,
using operators for notational simplicity, is

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

where ${\cal L}_{1}={\cal L}_{\rm cons}+{\cal L}_{\rm diff}+{\cal L}_{\rm
cap}+{\cal L}_{\rm ext}$ so that the operator ${\cal L}_{1}$ is the sum of a
conservative momentum transport term and diffusive, capillary and external
force terms. The first two terms are

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

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

${\bf D}=2\,\mu\,{{\bf S}\,,\quad{\bf S}=\tfrac{1}{2}\big{(}\nabla{\bf
u}+(\nabla{\bf u})^{T}\big{)}\,,}$ (7)

and $\mu$ is computed from $\chi$ using (1). The capillary term is

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

where $\sigma$ is the surface tension coefficient, $\nabla_{S}\sigma$ its
surface gradient, ${\bf n}$ the unit normal perpendicular to the interface,
$\kappa$ the sum of principal curvatures and $\delta_{S}$ a Dirac distribution
concentrated on the interface. Other equations are required to specify the
dependence of $\sigma$ from other physical quantities, such as temperature,
the presence of surfactants or electro-magnetic fields, therefore in the tests
section we will assume a constant $\sigma$. Finally ${\cal L}_{\rm ext}$
represents external forces. When the external force is gravity ${\cal L}_{\rm
ext}=\rho\,{\bf g}$.

The second, non conservative form, is

$\partial_{t}{\bf u}={\cal L}_{2}(\rho,{\bf u})-\frac{1}{\rho}{\nabla}p$ (9)

where ${\cal L}_{2}={\cal L}_{\rm adv}+({\cal L}_{\rm diff}+{\cal L}_{\rm
cap}+{\cal L}_{\rm ext})/\rho.$ The first term is

${\cal L}_{\rm adv}({\bf u})=-{\nabla}\cdot({\bf u}{\bf u})=-({\bf
u}\cdot\nabla){\bf u},$ (10)

and the other terms have been defined above.

The conservation of momentum, on an elementary control volume moving with the
interface $S$, leads to the following jump conditions across the interface

$\centering-\textlbrackdbl-p+2\mu\,{\bf n}\cdot{\bf S}\cdot{\bf
n}\textrbrackdbl_{S}=\sigma\kappa\@add@centering$ (11)

and

$\centering-\textlbrackdbl 2\mu\,{\bf t}^{(k)}\cdot{\bf S}\cdot{\bf
n}\textrbrackdbl_{S}={\bf t}^{(k)}\cdot\nabla_{S}\sigma\@add@centering$ (12)

where ${\bf t}^{(k)}$, $k=1,2$, are two independent tangent vectors and the
notation $\textlbrackdbl...\textrbrackdbl_{S}$ denotes the jump of a physical
quantity across the interface $S$. The Continuous Surface Force (CSF) [7]
formulation smoothens the surface tension force a few cells across the
interface and does not require the jump conditions (11,12), that are instead
implemented in the Ghost Fluid Method (GFM) [8].

###  2.2 Boundary conditions

A major difficulty with numerical simulations of fluid flow is the correct
implementation of the boundary conditions. In principle the conditions at
boundaries are well defined. For viscous, incompressible fluids we require
that the fluid sticks to the wall so that the fluid velocity there is equal to
the wall velocity

${\bf u}={\bf U}_{wall}.$

In a numerical setup, we can also impose periodic boundary conditions, as well
as inflow or outflow conditions (see TSZ).

###  2.3 Free-surface flow

Free-surface flow is a limiting case of flow with interfaces, in which the
treatment of one of the phases is simplified. For instance, for some cases of
air-water flow, we may consider the pressure $p$ in the air to depend only on
time and not on space (through, say, some function $p_{\rm air}(t)$) and the
viscous stresses in the air to be negligible. The jump conditions (11,12)
become boundary conditions on the border of the liquid domain

$\left.\left(-p+2\mu\,{\bf n}\cdot{{\bf S}}\cdot{\bf
n}\right)\right|_{S}={-}p_{\rm air}+\sigma\kappa$ (13)

and

$2\,\mu\,\left.{\bf t}^{(k)}\cdot{{\bf S}}\cdot{\bf n}\right|_{S}={\bf
t}^{(k)}\cdot\nabla_{S}\sigma\,.$ (14)

##  3 Numerical methods implemented in the code

###  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 cuboid grid, that can be defined as a cubic grid stretched
independently in the $x$, $y$ and $z$ directions, so that the centers of the
cells $\Omega_{i,j,k}$ are given by the intersection set of planes
$x=x_{i},y=y_{j},z=z_{k}$ and the cell boundaries are contained in the set of
planes $x=x_{i+1/2},y=y_{j+1/2},z=z_{k+1/2}$. However, in Paris the Front
Tracking method can indeed use stretched coordinates, while the routines
implementing the VOF method are still limited to cubic grids.

We use a finite volume discretization of the momentum equation and consider
staggered velocity and pressure grids. The staggered grid and control volumes
for the pressure $p$ are represented in Fig. 1. The corresponding control
volumes of the velocity component $u_{m}$ in direction $m$, $m=1,2,3$, are
shifted with respect to the control volume $\Omega_{i,j,k}$ surrounding the
pressure $p$. The control volumes for the velocity components $u_{1}$ and
$u_{2}$ in two dimensions, or for the corresponding momentum components, are
shown on Fig. 2. 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 TSZ for a more detailed discussion.) This
type of staggered representation is easily generalized to three dimensions.

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

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

| (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 indexes show the location of
the stored quantities. Thus, half indexes 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 denote the face of any control volume located in the positive or
negative Cartesian direction $m$, and ${\bf n}_{f}$ for the normal vector to
face $f$ pointing outwards of the control volume. On a cubic grid the spatial
step is $\Delta x=\Delta y=\Delta z=h$ and Eq. (15) becomes

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

where $u_{f}=u_{m\pm}={\bf u}\cdot{\bf n}_{f}$ is the velocity component
normal to face $f$.

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$,
$\mu$ and $C$.

###  3.2 Time Marching

Time marching can be performed in a first-order or second-order manner using a
small, possibly variable time step $\tau$ so that $t_{n+1}=t_{n}+\tau$. We
start with the description of the first-order time stepping. The interface is
first advanced in time as follows. In the Front-Tracking method, the front
points ${\bf x}_{k}$ are moved as

${\bf x}_{k}^{n+1}={\bf x}_{k}^{n}+\tau\,{\bf u}_{i}({\bf x}_{k}^{n},t_{n})$
(17)

where ${\bf u}_{i}({\bf x},t_{n})$ is an interpolation of the velocity field
at the location ${\bf x}$ (see Section 3.3 for the definition of the front
points). A marker function $I$ must then be constructed to compute the
material properties in the grid points, as described in Section 3.3.2.

In the VOF method the volume fraction field is updated as

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

where ${\cal L}_{VOF}$ represents the operator that updates the C data given
the velocity field ${\bf u}^{n}$ and is described in detail in Section 3.5.
Once the volume fraction $C$ (or the marker function $I$) is updated, the
material properties are computed from (1), while the interface geometry,
including its unit normal ${\bf n}$ and curvature $\kappa$, is estimated with
the different methods described in the next sections. The velocity field is
then updated with a predictor-corrector method. In the predictor step a
provisional velocity field ${\bf u}^{*}$ is computed. Two different versions
are used. One is the “non-consistent”, non conservative version of Eqs. (9,10)

$\displaystyle{\bf u}^{*}$ $\displaystyle=$ $\displaystyle{\bf
u}^{n}+\tau{\cal L}^{h}_{\rm adv}({\bf u}^{n})$ (19) $\displaystyle+$
$\displaystyle\frac{\tau}{\rho^{n+1}}\big{(}{\cal L}^{h}_{\rm
diff}(\mu^{n+1},{\bf u}^{n})+{\cal L}^{h}_{\rm cap}(C^{n+1})+{\cal L}^{h}_{\rm
ext}(C^{n+1})\big{)}\,,$

the other is a “mass-momentum consistent”, conservative version of Eqs. (5,6)

$\displaystyle\rho^{n+1}{\bf u}^{*}$ $\displaystyle=$
$\displaystyle\rho^{n}{\bf u}^{n}+\tau{\cal L}^{h}_{\rm cons}(\rho^{n},{\bf
u}^{n})$ (20) $\displaystyle+$ $\displaystyle\tau\big{(}{\cal L}^{h}_{\rm
diff}(\mu^{n+1},{\bf u}^{n})+{\cal L}^{h}_{\rm cap}(C^{n+1})+{\cal L}^{h}_{\rm
ext}(C^{n+1})\big{)}\,.$

Clearly the above operators depend on the discretization steps $\tau$ and $h$
as well as the fluid properties. The “mass-momentum consistent” advection
operator ${\cal L}^{h}_{\rm cons}$ is described in detail in [9]. In the
second step the pressure gradient corrects the velocity

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

This step constitutes the so-called projection method. The pressure is
determined by the requirement that the velocity at the end of the time step
must have zero divergence

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

which leads to an elliptic equation for the pressure

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

The whole set of operations above constitutes a first-order approximation in
time, which can be written

$({\bf f}^{n+1},{\bf u}^{n+1})={\bf L}({\bf f}^{n},{\bf u}^{n})$ (24)

where ${\bf f}^{n}$ is the interface data (either the front points ${\bf x}$
or the VOF fraction $C$) and the operator ${\bf L}$ consists in the steps
described above and is applied to the data at time $t_{n}$. A second-order
time scheme can be obtained by first computing a set of temporary variables
${\bf f}^{**}$ and ${\bf u}^{**}$

$({\bf f}^{**},{\bf u}^{**})={\bf L}({\bf f}^{n},{\bf u}^{n})$ (25)

followed by the update of the variables at time $t_{n+1}$ by the trapezoidal
rule

$({\bf f}^{n+1},{\bf u}^{n+1})=\frac{1}{2}\big{(}{\bf L}({\bf f}^{**},{\bf
u}^{**})+{\bf L}({\bf f}^{n},{\bf u}^{n})\big{)}$ (26)

where the operator ${\bf L}$ on the right hand side effectively is applied to
data at the intermediate time $t_{n+1/2}$. The Paris code implements both the
first-order and the second-order time schemes, controlled by the parameter
ITIME_SCHEME.

####  3.2.1 Non-conservative momentum advection

The non-conservative momentum advection in Eqs. (9, 10) amounts to integrate
over one time step the PDE

$\partial_{t}u_{m}={\cal L}_{{\rm adv},m}({\bf u})$ (27)

where

${\cal L}_{{\rm adv},m}={-}u_{j}\partial_{j}u_{m}$ (28)

for each value of the index $m$. Because of incompressibility (4) it is
equivalent to solve for a scalar field $\phi=u_{m}$ in the manner

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

Integrating the advection equation (29) in space, over a control volume
centered on a node of the scalar $\phi$, and in time one obtains

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

We use $F^{(\phi)}_{f}=\phi_{f}{\bf u}_{f}\cdot{\bf n}_{f}\tau/h$ as an
approximation of the flux on face $f$. Let $u_{f}={\bf u}_{f}\cdot{\bf
n}_{f}$. At first order the advecting velocity component $u_{f}$ is obtained
by simple averaging or centered interpolation. For example, we can consider
the first component of velocity $\phi=u_{1}=u$, whose control volume is
centered at $i+1/2,j,k$ (see Fig. 2). For the flux along the horizontal $x$
direction, the right face is located at index $i+1,j,k$. For the normal
component of the velocity on this face, the interpolated value is
$u_{f}=u_{1;i+1,j,k}=\frac{1}{2}(u_{1;i+1/2,j,k}+u_{1;i+3/2,j,k})$. On the
other hand for the flux along the $y$ direction, the top face is located at
index $i+1/2,j+1/2,k$ and the advecting velocity component is now
$u_{f}=u_{2;i+1/2,j+1/2,k}=\frac{1}{2}(u_{2;i,j+1/2,k}+u_{2;i+1,j+1/2,k})$.

Contrary to the estimates of $u_{f}$ above, the estimation of the advected
quantity $\phi_{f}$ may involve more complex and higher-order schemes. Indeed
one-dimensional interpolation schemes incorporating flux limiters are
typically used. Most of these schemes are described in TSZ together with their
usage in computing the face fluxes in the bulk. We describe their general
properties here shortly. Since the schemes are one-dimensional we can consider
a variable $\phi$ defined on a regular one-dimensional grid. Specializing
further the example, we assume that the variable $\phi$, akin to the velocity
component $u_{1}$, takes values $\phi_{i+1/2}$ at half-integer grid point
indexes. We need to estimate the flux at integer index points $x_{i}$, thus we
need to predict $\phi_{i}=\phi_{f}$. To this aim, an interpolation function is
defined that computes this value as a function of the value of $\phi$ at the
four nearest points, and in an upwind manner based on the sign of
$u_{i}=u_{f}$, that is given by the previously-defined centered interpolation.
The face value $\phi_{i}$ is then approximately given by

$\phi_{i}=f(\phi_{i-3/2},\phi_{i-1/2},\phi_{i+1/2},\phi_{i+3/2},u_{i})$ (31)

where the function $f$ is both of sufficiently high order and such as to limit
the flux. To express the function $f$, Paris offers a choice of the ENO,
QUICK, Superbee, WENO, first-order upwind, Verstappen, or BCG schemes.

When momentum is advected with the non conservative formulation (19), $\phi$
in Eq. (30) is taken equal to one of the velocity components $u_{m}$. This
method is available whether one uses the VOF method or the Front-Tracking
method. With the VOF method we have extensively used a combination of QUICK,
away from the interface, and a first-order upwind near the interface, and a
combination of Superbee slope limiter, away from the interface, and its
modified version near the interface. These combinations have been found
numerically more robust than others, and several simulation results are
presented in [9].

With the Front Tracking method we use typically ENO or QUICK, as they offer a
good combination of accuracy and stability.

####  3.2.2 Mass-momentum consistent momentum advection

When using VOF, another momentum advection method is available that is
consistent with VOF advection and which implements a conservative scheme of
the form (20). This means that the same advection method is used near the
interface for both the VOF volume fraction $C$ and the velocity ${\bf u}$. In
other words, when there is a density jump on the interface, the discontinuity
of the momentum density $\rho{\bf u}$ is advected exactly at the same velocity
as the discontinuity of the mass density $\rho$. This can be expressed by
saying that the momentum advection and the VOF advection are consistent. An
explicitly formulated criterion for consistency is the following: if the
velocity ${\bf u}$ is uniform, then $\rho{\bf u}$ remains exactly proportional
to $\rho$. This should happen even when $\rho$ is obtained from the VOF-
advection of $C$ using Eqs. (1) and (18) and $\rho{\bf u}$ is obtained from
the operator ${\cal L}^{h}_{\rm cons}$. Such a “VOF-consistent” method is used
since it has been empirically found by several authors that the non-consistent
advection of the previous section was often unstable at large density ratios,
while consistent methods are more stable [10, 11, 12, 13, 14, 15, 16, 17, 18,
19].

The Paris code implements a modification of the classical momentum-preserving
scheme proposed by [20] for the case of a staggered grid and Volume of Fluid
(VOF) method. The scheme needs to be modified from the one in the previous
section only near the interface. In particular, away from the interface, the
density is constant and the scheme in (30) is already conservative. Near the
interface, with $\phi=\rho u_{m}$, the momentum advection can be written, see
[9] for details, as

$\big{(}\rho u_{m}\big{)}^{n+1}_{i,j,k}-\big{(}\rho
u_{m}\big{)}^{n}_{i,j,k}=-\sum_{\rm{faces}\,f}u_{m,f}\,F^{(\rho)}_{f}\,+\,\sum_{l=1}^{3}\widetilde{u}_{m}\,C^{(\rho)}_{l}\,,$
(32)

where $u_{m,f}$ is the one-dimensional interpolation of the velocity component
$u_{m}$ on face $f$ and $F^{(\rho)}_{f}$ is the mass flux through the same
face. The last sum represents the compressional term and it is related to the
fact that each term in Eq. (16) can be different from zero even if the flow
has zero divergence. In particular this sum is equal to zero for the mass-
conserving VOF method of [21].

The one-dimensional interpolation schemes are different if they are applied
away from or near the interface. The combinations that are used with the mass-
momentum consistent advection have been already discussed in the previous
section.

####  3.2.3 Implicitation of the viscous terms

The operator ${\cal L}^{h}_{\rm diff}(\mu^{n+1},{\bf u}^{n})$ may be treated
in part implicitly. From (6,7) we have

${\cal L}_{\rm diff,j}(\mu,{\bf u})=({\nabla}\mu)\cdot\big{(}\nabla{\bf
u}+(\nabla{\bf u})^{T}\big{)}+\mu\nabla^{2}{\bf u}$ (33)

The first term on the RHS is left explicit, but the second term can be made
implicit by solving the linear problem

${\bf u}^{*}={\bf u}^{n}+\tau\,\mu^{n+1}\nabla^{2}{\bf u}^{*}\,.$ (34)

Then the discrete operator is defined as

${\cal L}_{\rm diff,j}^{h}(\mu^{n+1},{\bf
u}^{n})=({\nabla}\mu^{n+1})\cdot\big{(}\nabla{\bf u}^{n}+(\nabla{\bf
u}^{n})^{T}\big{)}+\big{(}{\bf u}^{*}-{\bf u}^{n}\big{)}\big{/}\tau$ (35)

where ${\bf u}^{*}$ is the solution of the linear problem (34). The
implicitation of the viscous terms is optional and controlled by a code
parameter.

###  3.3 Interface advection: Front-Tracking method

The interface advection in the “one-fluid” formulation is required to update
the material properties, see Eq. (1), and the interface geometry in order to
compute the surface tension forces. In the Paris code the interface can be
moved either by Front Tracking (FT) or by the VOF method. We describe the
former in this section. The FT method, in the context of simulations of two or
more immiscible fluids, refers to tracking the interface separating the
different fluids using moving connected marker points that represent the
interface. In our implementation the marker points are connected by triangular
elements, where the points are ordered in the same way for all elements,
allowing us to define an “inside” and an “outside” for each element. The
coordinates of the points are stored in arrays, in arbitrary order, with
separate integer arrays providing pointers to the previous and next points.
Thus, the points form a linked list where the location in the array provides
each point with a unique ID. The elements are stored in the same way, with
arrays containing pointers to the corner or node points. The coordinates of
the marker points are the main quantities stored for the points, but the
points also have arrays for various temporary quantities, such as velocities
and the surface force. The marker points and the elements connecting the
points together form the “front.” In addition to pointers to their corner
points, the elements also have pointers to the elements that share edges with
them. These are mostly used for modifications, or reconstruction, of the
front. Notice that generally only one front is needed, irrespectively of the
number of distinct interfaces, and that distinct interfaces, as present in a
simulation containing several bubbles and drops, can have different material
properties, for example surface tension. Figure 3 shows the layout of a
triangulated front separating two different fluids.

Figure 3: An interface separating two immiscible fluids represented by marker
points connected by triangular elements.

As the front is deformed, stretched and compressed by the flow, the size of
the elements changes as points move away from, or closer to, each other. We
keep the length of the edges of each element within about a quarter to a half
of the spacings of the fixed fluid grid, and to maintain that resolution,
points and elements are dynamically added and deleted. While many strategies
are possible, we add points by splitting the longest edge of an element by
adding one point and two elements, and delete points by collapsing the
shortest edge of an element, removing one point and two elements. The grid
quality can sometimes also be improved by changing the connectivity of the
elements but we generally find that doing so is not necessary. The addition
and deletion of front points and elements is shown in Fig. 4.

Topology changes in front tracking can be done in several slightly different
ways, but in all cases additional code is required to: (a) detect points where
coalescence/breakup should take place, and (b) restructure the front to
account for the topology change. For a short description of a topology change
algorithm and a few examples see [22] and for applications of a slightly
modified version of the algorithm to more complex flows see [23] and [24]. A
detailed description of a topology change algorithm is beyond the scope of the
current manuscript.

####  3.3.1 Connecting the front and the fluid grid

Since the Navier-Stokes equations are solved on a fixed grid, we have two
grids: the front/interface grid and the fixed grid. The motion of the
interface depends on the flow and the flow depends on where the interface is,
therefore information must be passed back and forth between the front and the
fixed fluid grid. To do so we need to identify what front point is close to
which fixed grid point and vice versa. For a regular structured grid, where
the grid lines are straight and evenly spaced, it is straightforward to locate
a point on the fixed grid that is closest to a given front point (using an INT
or a MOD function), but finding the front point closest to a given grid point
generally requires us to examine the distance to all the front points. Thus,
it is more efficient to do all communications between the front and the fixed
volume grid by looping over the front points. For periodic domains we allow
the front to move out of the domain resolved by the fixed fluid grid, and use
a MOD function to find the fixed grid point that would be closest to a given
front point if we moved the front back into the original domain.

Figure 4: Restructuring of a triangulated grid by adding and deleting points
and elements.

To transfer information between the fixed fluid grid and the moving front, we
need on one side to interpolate data from the grid to the front and on the
other side to spread, or “smooth,” data from the front to the fixed grid. The
first type of data transfer is required to move the front, where the velocity
at the front points is interpolated from the velocity on the fixed grid. On a
staggered grid each velocity component is interpolated separately. In general,
we have

$\phi^{l}_{f}=\sum_{ijk}w^{l}_{i,j,k}\,{\phi_{i,j,k}}$ (36)

Here, $\phi^{l}_{f}$ is a quantity, such as the velocity, on the front at
point $l$, $\phi_{i,j,k}$ is the same quantity on the fluid grid,
$w^{l}_{i,j,k}$ is the weight of each grid point with respect to front point
$l$ and the sum is over grid points “close” to the front points. Generally the
same time integration method is used for the advection of the points as is
used for updating the fluid velocities.

The second type of data transfer is usually referred to as smoothing, since it
replaces a quantity defined at a front point on the sharp interface, such as
surface tension, with a distribution on the fixed grid, where each fixed grid
point receives a value according to how close it is to the front point. Unlike
the interpolation of a quantity, such as the velocity from the fixed grid to a
front point, smoothing usually involves quantities like a force, that are
given in terms of force per unit interface area on the front but must be
converted to force per unit volume on the fixed grid, so that the total force
is conserved. Thus, the quantity smoothed must be scaled by the ratio of the
area associated with each front point divided by the volume of a fixed grid
cell

${\phi_{i,j,k}}=\sum_{l}\phi_{f}^{l}\,w_{i,j,k}^{l}\,{\Delta{A}_{l}\over{\Delta
x\Delta y\Delta z}},$ (37)

where $\Delta{A}$ is a surface area of a front element and $\Delta x$, $\Delta
y$ and $\Delta z$ are the grid spacings.

Several interpolation/smoothing functions can be used, but in the PARIS code
we use a smoother interpolation function originally introduced in [25] which
involves four grid points in each coordinate direction, or 64 points total.
The weights are given by

$w^{l}_{i,j,k}=d(r_{x})\,d(r_{y})\,d(r_{z}),$ (38)

where $r_{x}$ is the scaled distance (by the grid spacing) between $x^{l}_{f}$
and the grid line located at $x_{i}$; $r_{y}$ and $r_{z}$ are defined in the
same way. In our case, the weighting function $d(r)$ is given by

$d(r)=\left\\{\begin{array}[]{ll}(1/4)(1+\cos(\pi r/2)),&\left|r\right|<2,\\\
0,&\left|r\right|\geq 2.\end{array}\right.$ (39)

Using fewer points gives a sharper transition zone but sometimes leads to
wiggles, particularly for stiff problems. The interpolation function is
bounded with weights that sum to one in addition to having various desirable
symmetry properties. For a discussion see the reference above.

####  3.3.2 Constructing the marker function I

Once the front has been moved, a marker function must be constructed on the
fixed grid to assign the different material properties to each grid point.
This can be done in many different ways, but one of the consideration is that
fronts that are so close to each other that the flow between them is not
resolved, must be handled in a plausible way. Usually this means that the
marker function must retain its correct value on both sides of the double
front. In the PARIS code we do this by working with the gradient of the marker
function I, which in the limit of a sharp interface should be a delta function
defined only on the interface. The delta function is then treated in the same
way as the surface tension and smoothed onto the fixed grid. The grid value of
the components of the gradient is given by Eq. (37) of the previous section

$\big{(}{\nabla}I\big{)}_{i,j,k}=\sum_{l}\big{(}\Delta I\,{\bf
n}\big{)}_{l}\,w_{i,j,k}^{l}\,{\Delta{A}_{l}\over{\Delta x\Delta y\Delta
z}}\,,$ (40)

where $\Delta I$ is the jump in the value of the marker function across the
interface (usually a given constant value) and ${\bf n}$ the unit normal to
the element. Once the gradient field has been smoothed onto the fixed grid, we
can integrate it to recover the marker function I. For example, on a staggered
grid

$I_{i,j,k}=I_{i-1,j,k}+{\Bigl{(}{\partial I\over\partial
x}\Bigr{)}}_{i-1/2,j,k}\Delta x\,.$ (41)

To maintain symmetry, the marker at a given point is constructed as the
average of the integration from all the neighboring points, leading to a
linear system that is solved iteratively. Using standard second-order centered
finite-difference approximations, the linear system is an approximation to

$\nabla^{2}I=\nabla\cdot(\nabla I)_{f},$ (42)

where the subscript on the gradient of $I$ on the right hand side means that
it comes from the front. Since the right hand side is known, as it is deduced
from the position of the front through (41), this equation amounts to the
Poisson equation

$\nabla^{2}I=\nabla\cdot\mathbf{G}_{f}$ (43)

which must be solved to find $I$ for a given
$\mathbf{G}_{f}=\big{(}{\nabla}I\big{)}_{f}$. In the current version this
equation is solved for the entire grid to keep the implementation simple, but
this can easily be changed to involve only grid points next to the interface,
where the value of the marker function changes as the front moves. Integrating
the gradient of the marker function on the fixed grid is, of course, only one
way to construct it. We have, however, found that doing so generally leads to
a smooth but compact transition from one fluid to the other. In addition,
since the gradients of two interfaces bounding a very thin film cancel each
other when transferred to the fixed grid, the marker there will “disappear”.
This seems like the proper way to treat films too thin to be resolved on the
fixed grid.

Figure 5: Computation of the surface force on a triangulated grid by
integrating over the edges of an element.

###  3.4 Surface tension: Front-Tracking method

In simulations of flows with sharp interfaces the front serves two main
functions. The first is the advection of the marker function I, as described
above, and the second is the computation of the surface tension described
below. As with most of the other operations for the front, finding the surface
force can be done in several different ways. In the PARIS code we compute the
force on the front and transfer it to the fixed fluid grid, where it is added
to the discrete Navier-Stokes equations.

To ensure that the force is conserved as we transfer it from the front and
onto the fixed volume grid, we work with the total force on a small area. The
total force on a small region surrounding a front point is computed by
dividing each front element into three equal parts, each connected to one
nodal point, and computing the pull on their side as described below. The
force from each part is then added to the appropriate nodal point. When the
surface force on all the elements has been computed, it is transferred to the
fixed grid and converted into a force per unit volume by (37). Working with
the total force on a surface element, rather than the force per area, makes it
easier to ensure the total force is conserved when it is transferred between
the front and the fixed volume grid.

To find the surface force we use the fact that the total force on a surface
element can be found by integrating the “pull” on its edges

${\bf f}_{\sigma}=\int_{\Delta A}\sigma\,\kappa\,{\bf
n}_{e}\,dA=\oint_{C}\sigma\,{\bf p}\,dl$ (44)

Here the Stokes theorem has been used to convert the area integral over
$\Delta A$ into a line integral along its boundary $C$. The unit vector ${\bf
p}={\bf t}\times{\bf n}_{e}$ is tangent to the interface and perpendicular to
the boundary of the interface element (see Fig. 5). By keeping the surface
tension coefficient $\sigma$ under the integral sign we allow for variable
surface tension. The benefit of using this expression is that we only need to
approximate tangents on the surface, not the curvature, and that the pull on
the side of one surface element is equal and opposite to the pull on the
adjacent element. Thus, the surface force is conserved in the sense that an
integral over a surface patch consisting of several surface elements is
guaranteed to give the same result as an integral over the boundaries of the
whole patch. For constant surface tension coefficient, the integral over a
closed surface is, in particular, guaranteed to be zero.

Figure 5 shows schematically how the integration is done. We assume that the
elements are flat and that surface tension, given at the front points, can be
non constant. It therefore has to be interpolated when approximating the
integral. The force at point ${\bf x}_{0}$ is computed as the “pull” on the
edges of the gray patch, surrounding it. The contribution from the element
connecting points ${\bf x}_{0}$, ${\bf x}_{1}$ and ${\bf x}_{2}$ is found by
first splitting it in three parts by connecting the centroid ${\bf
x}_{c}=(1/3)({\bf x}_{0}+{\bf x}_{1}+{\bf x}_{2})$ to the midpoints of the
edges ${\bf x}_{01}=(1/2)({\bf x}_{0}+{\bf x}_{1})$ and ${\bf
x}_{02}=(1/2)({\bf x}_{0}+{\bf x}_{2})$ and then finding the force on those
edges by approximating the integral using a midpoint rule

$\Delta{\bf f}_{\sigma}\approx\sigma_{1}\Delta s_{1}{\bf
p}_{1}+\sigma_{2}\Delta s_{2}{\bf p}_{2}\,.$ (45)

Here, the pull on each element is the cross product of the outward unit
normal, ${\bf n}_{e}$, and the tangent vector to the edge, i.e. $\Delta
s_{1}{\bf p}_{1}={\bf n}_{e}\times({\bf x}_{01}-{\bf x}_{c})$ and $\Delta
s_{2}{\bf p}_{2}={\bf n}_{e}\times({\bf x}_{c}-{\bf x}_{02})$, where $\Delta
s_{1}$ and $\Delta s_{2}$ are the lengths of the edges, and $\sigma_{1}$ and
$\sigma_{2}$ are the surface tensions at the midpoint of the edges. The normal
${{\bf n}_{e}}$ is found by the normalized cross product of two of the tangent
vectors to the edges of the element. After straightforward algebra, we have

$\Delta{\bf f}_{\sigma}\approx\frac{{\bf
n}_{e}}{3}\times\Big{(}\sigma_{1}\big{(}{\bf x}_{2}-{1\over 2}({\bf
x}_{1}+{\bf x}_{0})\big{)}-\sigma_{2}\big{(}{\bf x}_{1}-{1\over 2}({\bf
x}_{2}+{\bf x}_{0})\big{)}\Big{)}$ (46)

where the interpolated surface tension at the midpoint of the edges is

$\sigma_{1}={1\over 2}\Big{(}\frac{1}{2}\big{(}\sigma({\bf x}_{0})+\sigma({\bf
x}_{1})\big{)}+\frac{1}{3}\big{(}\sigma({\bf x}_{0})+\sigma({\bf
x}_{1})+\sigma({\bf x}_{2})\big{)}\Big{)}$ $={1\over 12}\big{(}5\sigma({\bf
x}_{0})+5\sigma({\bf x}_{1})+2\sigma({\bf x}_{2})\big{)}$ (47)

and $\sigma_{2}$ is given by a similar expression. The forces from the other
elements connected to point ${\bf x}_{0}$ are found in the same way and added
to give the total force on the nodal point, that is then “smoothed” onto the
fixed grid.

###  3.5 Interface advection: VOF method

When the interface location is captured by the VOF method, a variable
$C_{i,j,k}$ is initialized. It is equal to the fraction of the cell
$\Omega_{i,j,k}$ that is filled with the reference fluid $1$. The
implementation of the VOF method is limited for the time being to cubic cells
of edge length $h$. Moreover, we use a rescaling of space and time variables
so that the cell size and the time step are both $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.5.1 Normal vector determination

The VOF method proceeds by a sequence of interface reconstructions and
advections of $C$. In the reconstruction step, one attempts to recover the
interface geometry from the VOF data $C_{i,j,k}$. In the Paris code we use
already-published methods (see for example TSZ) that have been experienced to
work satisfactorily. One first determines the interface normal vector ${\bf
n}$, then solves the problem of finding a plane perpendicular to ${\bf n}$
which cuts the cube with exactly the volume $C_{i,j,k}$. In Paris two methods
exist for normal vector determination. The most frequently used is the Mixed-
Youngs-Centered Scheme (MYCS) [26], described also in TSZ. However, when a
quick determination of the normal is needed but not very accurate, the finite
difference method ${\bf n}=\nabla^{h}C$ (the Youngs scheme [27]) can also be
used.

####  3.5.2 Plane constant determination

Once the the interface normal vector ${\bf n}$ is determined, a new, colinear
normal vector noted ${\bf m}$ and having unit “box” norm is deduced from ${\bf
n}$, that is $||{\bf m}||_{1}=|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$ (48)

cuts exactly a volume $V$ of the cube. The origin of the coordinate system is
taken at the corner of cubic cell $i,j,k$ with the smallest coordinate values.
The reader is reminded that we are using rescaled units of space, so that
$0\leq V\leq 1$ and $0\leq\alpha\leq 1$. Then $\alpha$ is determined by the
resolution of a cubic equation [28]. This resolution, and similar ones often
used in the VOF method and derived in [28], are implemented in a kind of small
library contained in the single file vof_functions.f90.

####  3.5.3 Volume initialization

Before any VOF interface tracking is performed, the field of $C_{i,j,k}$
values must be initialized. The Paris code avoids inaccurate initializations
that for example initialize a sphere as a set of $C_{i,j,k}$ values which are
either $0$ or $1$, a so called “staircase” or “lego” initialization. There are
two ways in which initialization can be improved over the lego one. In the
“subgrid” initialisation, the mesh cells are subdivided into $n_{I}^{3}$
subcells (where $n_{I}$ is a tunable parameter, called REFINEMENT in the
code). Then a “lego” initialization is performed trivially in the subcells.
For example, if the initial interface is defined implicitly by the equation
$\phi({\bf x})=0$ where $\phi$ is a smooth implicit function (akin to a level-
set function) then the trivial “lego” initialization in each subcell is
$c_{I}=\chi\big{(}\phi({\bf x}_{c})\big{)}$, where $\chi$ is the Heaviside
function and ${\bf x}_{c}$ the subcell center. The cell value $C_{i,j,k}$ is
then given by the sum of the $c_{I}$ values divided by the subcells number
$n_{I}^{3}$. In tests it was found that $n_{I}=8$ was sufficient. However
$n_{I}\geq 8$ leads to a very large number of evaluations of the function
$\phi$ and a slow initialization. In order to avoid this, the Paris code may
be linked to the Vofi library described in [29] and [30]. Then the
initialization is performed using highly-accurate numerical integration to
compute the fluid volumes in each grid cell, implicitly defined by the
equation $\phi({\bf x})<0$. The code is linked to the Vofi library with the
shell variable HAVE_VOFI set before compilation.

####  3.5.4 General split-direction advection

The reconstruction at time $t_{n}$ provides the approximate position of the
interface, which is used to compute the reference phase fluxes across the cell
boundary in order to update the volume fractions $C_{i,j,k}$ at time
$t_{n+1}$. The Paris code contains two methods for split-direction advection,
which can be selected by the user: Lagrangian Explicit (LE) advection, with
the keyword VOF_ADVECT = LE, and Weymouth and Yue (WY) advection, with the
keyword VOF_ADVECT = WY. The LE advection is also called “Calcul d’Interface
Affine par Morceaux” (CIAM) which is French for “Piecewise Linear Interface
Calculation” (PLIC), however PLIC refers to generic VOF methods with a
piecewise linear reconstruction step, while CIAM refers to a specific type of
advection method first described in the archival literature in [31] and
classified as the “LE” method in [32]. The main advantage of both LE and WY is
that they avoid overshoots ($C_{i,j,k}>1$) and undershoots ($C_{i,j,k}<0$).
Moreover WY conserves mass to machine accuracy. These methods are described in
detail in TSZ for LE/CIAM and in [21] for WY. The reader may also refer to [9]
for a condensed description of both methods.

An important operation in some simulations is the “clipping” procedure, and a
small parameter $\epsilon_{c}$ is defined to this purpose. After advection,
all cells that have $C_{i,j,k}<\epsilon_{c}$ are set to 0 and all cells that
have $C_{i,j,k}>1-\epsilon_{c}$ are set to 1. This removes some, but not all,
of the wisps, floatsam and jetsam that are generated. In the current version
of the code the default value is $\epsilon_{c}=10^{-8}$. This is a rather high
value used mainly with WY in simulations at very high density contrast. The
code has been observed to run well also with the smaller threshold value
$\epsilon_{c}=10^{-12}$ when the ratio of physical parameters is less extreme,
or with no threshold at all with the LE advection. Other VOF schemes, and in
particular unsplit schemes, are not reported to require clipping.

###  3.6 Surface tension: VOF method

####  3.6.1 CSF method

For simplicity, we consider only the case where $\sigma$ is constant. In the
Continuous Surface Force (CSF) method [7] the force $\sigma\kappa{\bf
n}\delta_{S}$ in the capillary term (8) is written as

$-\sigma\kappa\nabla\chi=-\sigma\kappa^{h}\nabla^{h}C.$ (49)

where we have used the properties of the Heaviside function $\chi$. One of the
advantages of this formulation is that it is a “well-balanced” method (see
TSZ, Chapter 7 or reference [33]).

An approximation for $\kappa$ needs to be found to use the CSF method. A good
estimate is obtained using so-called height functions.

####  3.6.2 Height Functions

We give some details about height functions since it is a relatively novel
aspect of the code. Height functions were introduced in [34] and further
discussed, tested and improved in several papers [35, 36]. A height function
is a function on the discrete grid that gives the elevation of the interface
with respect to a reference plane. The use of height functions greatly
improves the accuracy of VOF methods since it allows us to neglect small
inconsistencies in the VOF representation. On one hand a small VOF floatsam in
a cell has only a very small influence on the height-function calculation. On
the other hand if taken as an indication of the presence of the interface in
the given cell it would create a large error on the interface location (see
Fig. 6a). With high resolution, height functions can also yield the position
of the interface to fourth-order accuracy [35].

Figure 6: (a) A small floatsam $F$ in a cell away from the interface is
negligible when height functions are used to determine the location of the
interface. On the other hand requiring the interface to pass near each
centroid $C_{i}$ has a large effect. A small inconsistency such as near point
$I$ is also ignored by the height function. (b) Two cases where the height
function expression (50) is appropriate with four cells ($n_{c}=3$, see text).
For more vertical lines or larger curvatures the line exits the $4\times 1$
stencil through the top and bottom and the HF method cannot be used. (c) To
check the validity of a the HF calculation one needs to have one full cell
($C=1$) below and one empty cell ($C=0$) above, or the converse.

We consider for simplicity the case of an approximately horizontal interface,
then the reference plane is aligned with the $x,y$ plane. A local height $H$,
rescaled by the grid size $h$, may be defined as

$H=\sum_{p=k_{0}}^{p=k_{0}+n_{c}}C_{i,j,p}.$ (50)

and is illustrated on Fig. 6b. The expression (50) defines a “vertical” height
with the reference plane passing through point $O$ and the base of the bottom
cell on Fig. 6b. More general cases, with arbitrary orientation of the
interface, are considered in A.

The computation of the curvature may then be performed by reconstructing a
polynomial approximation to the height function. To illustrate this, we take
again the case of an approximately horizontal interface. Then we fit the local
heights by the function

$H(x,y)=\frac{a_{1}}{2}x^{2}+\frac{a_{2}}{2}y^{2}+a_{3}xy+a_{4}x+a_{5}y+a_{6}$
(51)

where the coefficients $a_{i}$ are computed from the heights data using finite
differences (see Eq. (80)). The curvature is then

$\kappa=\epsilon\frac{a_{1}(1+a_{5}^{2})+a_{2}(1+a_{4}^{2})-2a_{3}a_{4}a_{5}}{(1+a_{4}^{2}+a_{5}^{2})^{3/2}}$
(52)

where $\epsilon=1$ if the interface is in the “canonical” position (normal
pointing upwards). The above method is possible only if all nine vertical
heights are available in the $x,y$ plane. If they are not, fallback methods
are used, details of which are given in A.

The performance of our method for the computation of curvature is shown in
Figs. 7 and 8. The error was computed for a collection of diameter-to-cell-
size ratios $D/h$. For each value of $D/h$ the error computation was repeated
for an ensemble of $N$ sphere centers located randomly. The $L_{\infty}$ norm
for a given sphere center is the maximum difference between the sphere
curvature and the numerically obtained curvature. The error reported on these
two figures is the maximum $L_{\infty}$ error for the whole ensemble of $N$
spheres. We checked that the error varies little when $N$ is increased above
16, and the standard test case for curvature distributed with the code uses
this value of $N$. The differences between the error norms obtained for
Basilisk and Paris are discussed in A.

Figure 7: Maximum $L_{\infty}$ error norm in two dimensions for the curvature
estimated for a cylinder using the height function method in Paris simulator
and Basilisk. The mixed-height option is set in both codes. Figure 8: Maximum
$L_{\infty}$ error norm in three dimensions for the curvature estimated for a
sphere using the height function method in Paris and Basilisk. The mixed-
height option is set in both codes.

###  3.7 Pressure solver

####  3.7.1 In-code Gauss-Seidel solver

The default Poisson solver used to invert the elliptic operators appearing in
Eqs. (23) and (34) is a red-black Gauss-Seidel (GS) solver with overrelaxation
[37]. Both equations can be discretized in the form

$\displaystyle A_{1,i,j,k}\,p_{i-1,j,k}$ $\displaystyle+$ $\displaystyle
A_{2,i,j,k}\,p_{i+1,j,k}+A_{3,i,j,k}\,p_{i,j-1,k}+$ $\displaystyle
A_{4,i,j,k}\,p_{i,j+1,k}$ $\displaystyle+$ $\displaystyle
A_{5,i,j,k}\,p_{i,j,k-1}+A_{6,i,j,k}\,p_{i,j,k+1}-$ $\displaystyle
A_{7,i,j,k}\,p_{i,j,k}$ $\displaystyle=$ $\displaystyle A_{8,i,j,k}\,,$ (53)

where the coefficients verify the relations

$\displaystyle A_{7,i,j,k}$ $\displaystyle=$
$\displaystyle\sum_{p=1}^{6}A_{p,i,j,k}$ (54) $\displaystyle A_{2,i,j,k}$
$\displaystyle=$ $\displaystyle A_{1,i+1,j,k}$ (55) $\displaystyle
A_{4,i,j,k}$ $\displaystyle=$ $\displaystyle A_{3,i,j+1,k}$ (56)
$\displaystyle A_{6,i,j,k}$ $\displaystyle=$ $\displaystyle A_{5,i,j,k+1}$
(57)

and are constructed by interpolations of $1/\rho$, for Eq. (23), or $\mu$, for
the implicit part (34) of the viscous term in the momentum equation. The GS
solver iterates the assignment

$\displaystyle p_{i,j,k}$ $\displaystyle\leftarrow$
$\displaystyle\left(1-\beta\right)p_{i,j,k}+\frac{\beta}{A_{7,i,j,k}}\Big{(}A_{1,i,j,k}\,p_{i-1,j,k}$
(59) $\displaystyle+$ $\displaystyle
A_{2,i,j,k}\,p_{i+1,j,k}+A_{3,i,j,k}\,p_{i,j-1,k}+A_{4,i,j,k}\,p_{i,j+1,k}$
$\displaystyle+$ $\displaystyle
A_{5,i,j,k}\,p_{i,j,k+1}+A_{6,i,j,k}\,p_{i,j,k-1}-A_{8,i,j,k}\Big{)}$

where $\beta$ is an overrelaxation parameter that can be set by the user. A
value of $\beta=1.3$ is typically used.

####  3.7.2  library multigrid solver

The Hypre library, developed at the Lawrence Livermore National Laboratory
(LLNL), is also an option to solve the elliptic equations with multigrid
iterative methods. Since a structured grid is used in Paris, a few solvers of
the Hypre library can be used. The SMG and PFMG multigrid solvers have been
implemented in the code and used for large-scale simulations using up to
64,536 cores. Both SMG and PFMG are parallel semicoarsening multigrid solvers.
The difference lies in that the SMG solver uses plane smoothing while the PFMG
solver uses pointwise smoothing. The plane-smoothing feature makes the SMG
solver more robust but less efficient. In fact, the scaling performance of the
PFMG solver is much better than the SMG solver, since the smoothing operations
only involve a local stencil.

In order to take advantage of the higher efficiency of PFMG and the robustness
of SMG, a solution strategy has been implemented in the code. The PFMG solver
is used by default, however if the iteration diverges or fails to converge
within the maximum iteration number, the code will switch to the SMG solver
and redo the iteration. If the iteration converges, then the code will switch
back to PFMG for the next time step. For a large-scale simulation that runs
for a long time, this strategy has been shown to achieve a good performance,
balancing robustness and efficiency.

The Hypre library, at least in the versions we use, appears to control the
tolerance on the residual using the $L_{1}$ norm. The code then recomputes the
residual norms and controls the accuracy using a norm of the residual chosen
by the user among $L_{1}$, $L_{2}$ and $L_{\infty}$.

####  3.7.3 In-code multigrid solver

The code has also a native implementation of a multigrid solver for structured
grids with $2^{n}$ number of points per coordinate direction. In particular
the V-Cycle scheme is implemented and fully parallelized [37]. Relaxation
operations are applied starting first from the finest to the coarsest grid,
and then from the coarsest to the finest one, the number of relaxation
operations being a user-adjustable parameter. One advantage of having a native
multigrid solver is that it allows for an efficient solution of the Poisson
equation without the necessity of having external libraries installed in the
system, such as Hypre. Especially when running heavy three-dimensional
simulations in parallel the use of this native solver has been shown
advantageous in some systems with respect to Hypre in terms of memory
manipulation.

####  3.7.4 GPU-accelerated solver

A GPU-based solver is also available for solving the Poisson equation when a
significant number of iterations is required to achieve convergence. The
pressure is solved using a Jacobi method for Eq. (59). The need for the Jacobi
method instead of a Gauss-Seidel arises because of the intrinsic nature of GPU
devices. The usual domain decomposition parallelization allows the
implementation of the iterative step by using a simple for loop over the
indexes $i,j,k$. The sequentiality of the indexes cannot be achieved on GPU
devices, as in this case each index combination is ideally computed
simultaneously. In this sense, the larger number of iterations required by the
Jacobi method is mitigated by the speed-up provided by GPUs.

The memory handling is a critical aspect in GPUs applications and it is even
more critical in DNS. Although a Jacobi method intrinsically requires doubling
the memory usage for the pressure matrix $p$, it also enables the leanest data
transfer between CPUs and GPUs, which is a critical aspect of normal CUDA
applications. A Gauss-Seidel red-black solver is in principle possible and
could be beneficial for certain applications. In fact, let us assume that the
size of the matrix $p$ is $N_{p}=n_{x}n_{y}n_{z}$, the normal implementation
of a red-black Gauss-Seidel solver would be

for all $\Omega_{i,j,k}$ cells of ”red” type do

compute $p_{i,j,k}$ using (59)

end for

for all $\Omega_{i,j,k}$ cells of ”black” type do

compute $p_{i,j,k}$ using (59)

end for

check convergence

which inherently reduces the memory usage required. On the other hand, the
first for loop is not parallelizable in an efficient way in CUDA. Therefore,
such an algorithm would be beneficial from a memory standpoint, but would
improve the computational time only if $N_{p}$ is at least 4 times greater
than the number of GPU process available. As this is usually not the case, the
beneficial effects of a red-black algorithm are limited, although it will be
object of future studies.

The implementation of the algorithm is achieved by means of the open-source
CUDA library for C developed by NVIDIA, while the intercommunication between
processors is still achieved by using MPI. For this reason, an interface
between Fortran90 and C is created in module_CUDA.f90. By passing through the
interface, each process transfers the coefficients matrix $A$ and the pressure
matrix $p$ to the C/CUDA environment (in poissonCUDA.cu) where the iteration
step is performed. The boundary conditions, as well as the MPI communication,
are enforced in the environment that originally created the MPI communication,
hence these functions are programmed in cudaFun.f90.

####  3.7.5 Free-surface flow solver

A free-surface flow solver is implemented in Paris, which is designed to apply
a free-surface condition as described in Section 2.3. The implementation of
the boundary conditions (13,14) is here limited to inviscid flows, hence only
the jump condition (13) along the normal direction is considered. The free-
surface solver uses the VOF method implemented in Paris to track the
interface. The flow is then solved only in the phase with $\chi=0$, using the
same numerical methods described previously. For the purpose of description
the solved phase will be called liquid and the unsolved phase will be gas. In
the region occupied by the gas phase the conservation equations are not
solved, instead a homogeneous pressure field is assumed, that can vary only in
time. The time variation of the pressure is determined by a polytropic gas law

$p_{c}=p_{0}\left(\frac{V_{0}}{V_{c}}\right)^{\gamma}\,,$ (60)

where $V_{c}$ is the total volume of the gas pocket at pressure $p_{c}$,
$p_{0}$ and $V_{0}$ are respectively the reference pressure and volume of the
gas phase, and finally $\gamma$ is the heat capacity ratio. The gas phase
pressure along with the pressure jump due to surface tension is then applied
as a Dirichlet boundary condition for the pressure in the liquid flow, as
given in Eq. (13).

The method used to apply this pressure boundary condition is inspired by the
idea of Fedkiw and Kang [8, 38], often referred to as the Ghost Fluid Method
(GFM). Let the time-varying pressure in the unsolved phase be $p_{c}$. Special
care is required in the discretization of the elliptic equation (23) for
liquid cells near the interface. Cells that contain mostly gas are excluded
from the solution, so that only cells where $C<0.5$ are solved (we recall that
the convention is to have $C=0$ in the liquid). Figure 9 shows a
representation of a 2D grid with a section of an interface. The grey area
represents a gas-filled volume. Cells that contain a filled circle are
included in the pressure solution, while cells without a marker in the center
are excluded.

Figure 9: A 2D section of the numerical grid, showing part of a gas bubble in
grey. Circles represent computational cell nodes where pressure is calculated.
Triangles indicate scalar velocity components on the computational cell faces.
Filled triangles indicate values which are found by solving the governing
equations, while unfilled triangles represent boundary values found by
extrapolation.

Since only the liquid phase is solved, only the liquid density is applied.
Furthermore, for cubic cells $\Delta x=\Delta y=\Delta z=h$, where $h$ is the
constant grid spacing.

$p_{i,j}$

$p_{i-1,j}$

$p_{s,i,j-1}$

$p_{i,j+1}$

$\delta_{i,j-1/2}$

$\delta_{i+1/2,j}$

$p_{s,i+1,j}$

Figure 10: Discretisation of the pressure equation near the interface

The stencil for the pressure gradient components has to be changed near the
interface when a neighbouring cell falls inside the gas phase. The pressure in
this cell must be substituted by a surface pressure $p_{s}$. We apply a
finite-difference gradient approach as Chan [39]. As an example, the
approximation for the pressure gradient components for the cell with indexes
$i$ and $j$ in Fig. 10 is written

$\nabla^{h}_{x}\,p_{i+{1}/{2},j}=\dfrac{p_{s,i+1,j}-p_{i,j}}{\delta_{i+{1}/{2},j}}\,\,;\quad{\nabla}^{h}_{y}\,p_{i,j-{1}/{2}}=\dfrac{p_{i,j}-p_{s,i,j-1}}{\delta_{i,j-{1}/{2}}}\,,$
(61)

where $\delta$ is the distance between the pressure node under consideration
and the intersection with the interface. The surface pressure $p_{s}$ on the
liquid side of the interface is found by adding to the gas pressure $p_{c}$
the Laplace pressure jump. The pressure $p_{c}$ inside each cavity is known
from the polytropic law (60). The interface pressure $p_{s}$ in the first
expression of (61) will then be

$p_{s,i+1,j}=p_{c,i+1,j}+\sigma\frac{\kappa_{i,j}+\kappa_{i+1,j}}{2}\,.$ (62)

From Eqs. (61) and (62) it is clear that an accurate computation of the
interface curvature as well as an accurate prediction of the interface
location are important parameters to ensure the accuracy of the pressure
solution. Since the height function is the approximate interface distance from
some reference cell boundary in a given direction, it is used for $\delta$.
When the interface configuration is such that a height cannot be obtained in
the required direction, the distance is approximated by using a plane
reconstruction of the interface in the staggered volume.

##### Extrapolation of the velocity field

The resolved velocity components right next to the interface will require
neighbours in the gas phase to discretize the momentum advection term. These
values in the gas phase can be seen as boundary values to the resolved
velocities. In order to find neighbours in the gas phase, we extrapolate the
resolved velocities similarly to Popinet [40].

After calculating the liquid velocities using standard methods in Paris, the
boundary velocities in the gas phase are updated for the next time step from
the closest two velocity neighbours using a linear least-square fit. Let us
assume that the velocity field can be described as a linear combination

${\bf u}\left({\bf x}\right)={\bf A}\cdot({\bf x}-{\bf x}_{0})+{\bf u}_{0}$
(63)

where the components of the tensor ${\bf A}$ and of the vector ${\bf u}_{0}$
are the unknowns. In two dimensions we take a $5\times 5$ stencil around the
unknown gas velocity at location ${\bf x}_{0}$ and find the extrapolated
velocity ${\bf u}_{0}$ by minimizing the functional

$\mathcal{L}=\sum_{k=1}^{N}\big{|}{\bf A}\cdot({\bf x}-{\bf x}_{0})+{\bf
u}_{0}-{\bf u}_{k}\big{|}^{2}$ (64)

This is done first for all locations closest to the resolved velocities ${\bf
u}_{k}$ (“first neighbours”), whereafter the process is repeated for the
“second neighbours”. Note that only resolved velocity components are included
in the cost function, therefore the number of the known velocities $N$ can
vary depending on the shape of the interface. Furthermore, because of the
staggered grid, only one component of the velocity ${\bf u}_{0}$ is computed
at any location ${\bf x}_{0}$.

##### Ensuring volume conservation

An additional step is required to ensure that the extrapolated velocities are
divergence-free. This is required to ensure that the advection of $C$ is
conservative.

A similar approach to Sussman [41] is used. Only the first two layers of cells
inside the gas phase are considered and all other cells are disregarded. A 2D
example is presented in Fig. 11.

Liquid

Gas

Figure 11: 2D example of the problem to correct the extrapolated velocities
(unfilled triangles). A Poisson’s problem is solved in the cells marked with
an unfilled circle.

Similar to the projection step explained earlier, a “phantom” pressure is
obtained in these cells by solving a Poisson’s equation

$\nabla^{h}{\cdot}\bigg{(}\nabla^{h}\hat{P}\bigg{)}=\nabla^{h}{\cdot}{\tilde{{\bf
u}}}\,,$ (65)

where $\hat{P}$ is the “phantom” pressure and ${\tilde{{\bf u}}}$ is the
extrapolated velocity on the faces of the first two gas neighbours. $\hat{P}$
is only calculated in the cells represented by unfilled circles in Fig. 11. On
the liquid side of these cells, the solved velocities (filled triangles) are
used as a velocity boundary condition with the pressure gradient on this face
set to zero. On the gas side, a fixed pressure is prescribed in the cells (red
filled circles) outside the region where the “phantom” pressure is computed.
Only the extrapolated velocities (unfilled triangles) are then corrected by
the solved pressure gradient $\nabla\hat{P}$

${\tilde{{\bf u}}}^{n+1}={\tilde{{\bf u}}}-\nabla^{h}\hat{P}$ (66)

to ensure a divergence-free velocity field in the first two layers of cells
just inside the gas.

For more details on the numerical method and its application in idealized
microspalling, see [42] and [43]. This process is found in metals that under
shock loading melt and where, with the reflection of the shock wave from the
material free surface, cavities can nucleate, grow and merge.

###  3.8 Solid boundaries

Solids are defined in a “block” or “Lego” manner. A domain-wide binary flag
$s_{i,j,k}$ is defined that is equal to 0 inside the solid and 1 outside. A
set of link-based flags $s^{x}$, $s^{y}$ and $s^{z}$ is also defined (a link-
based array is data located on the velocity component collocation points, such
as $i+1/2,j,k$ for the horizontal velocity component $u$). The following
pseudo code is executed at initialisation relatively to the $x$ direction,
with $i=1,2,...,n_{x}$

for all $j,k$ do

$s^{x}_{1/2,j,k}\leftarrow s_{1,j,k}$

end for

for all $i,j,k$ do

$s^{x}_{i+1/2,j,k}\leftarrow s_{i,j,k}$

end for

and similarly for the $y$ and $z$ directions. The indexes $s^{x}$, $s^{y}$ and
$s^{z}$ are then used to “block” the velocity and the pressure correction on
the solid region and its boundary. This is done each time the velocity is
updated

for all $i,j,k$ do

$u_{i+1/2,j,k}\leftarrow{s^{x}_{i+1/2,j,k}}\,u_{i+1/2,j,k}$

end for

and similarly for the velocity components $v$ and $w$. The pressure correction
should not change the velocity on the solid boundaries, so on the links a
Neuman boundary condition for the pressure is established, which is equivalent
to setting to zero some coefficients

for all $i,j,k$ do

$A_{1,i,j,k}\leftarrow s^{x}_{i-1/2,j,k}\,A_{1,i,j,k}$

$A_{2,i,j,k}\leftarrow s^{x}_{i+1/2,j,k}\,A_{2,i,j,k}$

(similarly for $A_{3},A_{4}$ and $s^{y}$, and for $A_{5},A_{6}$ and $s^{z}$)

$A_{7,i,j,k}\leftarrow\sum_{p=1}^{6}A_{p,i,j,k}$

end for

The solid domain can be initialized by implicit functions or by loading a file
containing the $s_{i,j,k}$ data. In both cases it is important that the
presence of the solid does not make the linear system (23) ill-posed. This
will happen for example if the boundary conditions are Dirichlet for the
velocity at the entry of a channel and the solid completely blocks the
channel. To avoid this type of problem, there is a small utility program
“rockread.c”, distributed with the code, that checks that the fluid
“percolates”.

###  3.9 Lagrangian point-particle model

A Lagrangian point-particle (LPP) model has been implemented in PARIS, that is
fully coupled to the VOF method to provide a multi-scale modeling strategy to
simulate liquid atomization: the large-scale interfaces on the bulk liquid are
resolved by VOF, while the small droplets are represented by the LPP model.
The details of the model can be found in [44] and only a brief summary is
given here.

####  3.9.1 Overview

The combined VOF-LPP algorithm is shown in Fig. 12. The Navier-Stokes
equations and the advection equation for the volume fraction $C$ are solved
for the resolved flow. Then the $C$ field is tagged to identify different
liquid structures. The tagging approach of [45] is used, and cells that are
attached to each other, with respect to the liquid phase, will have the same
tag number. The size, aspect ratio, and distance from the bulk flow of each
liquid structure are computed in the tagging process. Small resolved droplets
(RD) that satisfy the LPP criteria will be submitted to the RD-to-LPP
conversion routine. The motion equation is then solved for each LPP droplet to
update its velocity and position. Afterwards, the droplet is examined to
decide whether or not it should be converted back to a resolved droplet.

LPP droplets and the resolved flow (RF) are two-way coupled. The RF
properties, such as fluid velocity, are needed to calculate the forces acting
on each LPP droplet. On the other hand, these forces need to be subtracted
from the momentum equation of the flow, thus appearing as an additional source
term.

Figure 12: Flow chart of the combined VOF-LPP algorithm.

####  3.9.2 Two-way coupling between LPP and RF

The non-conservative momentum equation (9,10) is now written as

$\rho\big{(}\partial_{t}{\bf u}\,+\,{\bf u}\cdot\nabla{\bf
u}\big{)}=-{\nabla}p+{\nabla}\cdot{\bf D}+{\sigma\,\kappa\,\delta_{S}\,{\bf
n}}-{\bf f}_{p}$ (67)

where ${\bf f}_{p}$ is the closure term that accounts for the backward effect
of the LPP droplets on the resolved flow. The LPP model approximates each
small droplet as a point mass, hence the droplet-scale flow is not resolved.
To accurately track the LPP droplets in the Lagrangian framework, the force
exerted on each droplet is calculated in terms of the undisturbed flow field
properties. The closure is typically given by the force model or the so-called
equation of motion (EOM) [46]

$\displaystyle\frac{d{\bf x}_{p}}{dt}$ $\displaystyle={\bf u}_{p}\,,$ (68)
$\displaystyle\frac{d{\bf u}_{p}}{dt}$ $\displaystyle=\frac{{{\bf u}}-{\bf
u}_{p}}{\tau_{p}}\phi+\frac{\rho}{\rho_{p}}\frac{D{{\bf
u}}}{Dt}+\frac{C_{m}\rho}{\rho_{p}}\left(\frac{D{{\bf u}}}{Dt}-\frac{d{\bf
u}_{p}}{dt}\right)+{\bf g}\,,$ (69)

where $\rho$ and ${{\bf u}}$ are the density and velocity of the undisturbed
flow, and ${\bf x}_{p}$, ${\bf u}_{p}$,
$\tau_{p}={\rho_{p}\,d_{p}^{2}}/(18\,\mu)$, $\rho_{p}$, and $C_{m}$ the
position, velocity, response time, density, and added-mass coefficient of the
LPP droplet. A tri-linear interpolation is used to interpolate the fluid
velocity from the grid to the position of the LPP droplets. The modified
gravity acceleration is denoted by ${\bf g}=(1-\rho/\rho_{p}){\bf
g}^{\prime}$, where ${\bf g}^{\prime}$ is the gravity acceleration. For the
time integration of Eqs. (68) and (69) a second-order predictor-corrector
method is used, which is consistent with the algorithms used for the resolved
flow.

The force terms on the right hand side of Eq. (69) represent the quasi-steady
force, the pressure-gradient force, the added-mass force, and the gravity
force, respectively. The Basset-history force has been ignored for simplicity,
while the Faxén and lift forces have been neglected because the LPP model is
here applied to small droplets, and the effect of inhomogeneous ambient flows
on the droplet force is expected to be small. The quasi-steady force is
expressed as the Stokes drag multiplied by a correction term $\phi$, which is
a function of the Reynolds number $\mathrm{Re}_{p}={\rho\,d_{p}\,|{{\bf
u}}-{\bf u}_{p}|}/{\mu}$ [47],

$\displaystyle\phi(\mathrm{Re}_{p})$
$\displaystyle=1+0.15\,\mathrm{Re}_{p}^{0.687}+0.0175\,\mathrm{Re}_{p}\left(1+\frac{42500}{\mathrm{Re}_{p}^{1.16}}\right)^{-1}\,.$
(70)

In atomization, the inter-droplet interaction can be ignored, as the small
droplets are quickly dispersed away from the bulk liquid and the related
volume fraction is relatively small.

As a consequence of Newton’s third law, the force ${\bf f}_{p}$ exerted on the
LPP droplets needs to be subtracted from the resolved flow and is referred to
as backward coupling. This force in the momentum equation (67) is calculated
as

$\displaystyle{\bf f}_{p}=\sum_{i=1}^{N_{p}}\mathbf{F}_{fp,i}\,G({\bf x}-{\bf
x}_{p,i})\,,$ (71)

where $N_{p}$ is the total number of LPP droplets. The Gaussian function
$G({\bf x}-{\bf x}_{p,i})$ is a numerical representation of the LPP droplet
coupling force [48]

$\displaystyle G({\bf x})=(2\pi\mathcal{L}^{2})^{-3/2}\exp(-|{\bf
x}|^{2}/2\mathcal{L}^{2})\,,$ (72)

where $\mathcal{L}$ controls the size of the region where the force should be
distributed. Note that only the force due to fluid-LPP-droplet interaction
needs to be subtracted from the momentum equation

$\displaystyle\mathbf{F}_{fp,i}=m_{p,i}\left(\frac{d{\bf u}_{p,i}}{dt}-{\bf
g}\right)\,.$ (73)

####  3.9.3 Two-way conversion between LPP and resolved droplets

The droplets that are generated in the atomization process are converted to
LPP droplets when they satisfy the following criteria.

##### LPP Criteria

The criteria to determine whether a resolved droplet (RD) should be
represented as a LPP droplet are based on its volume $V_{p}$, aspect ratio
$\gamma_{p}$, and position ${\bf x}_{p}$. $V_{p}$ is required to be smaller
than a critical value $V_{crit}\simeq(4\,h)^{3}$, and $\gamma_{p}$ close to
one. Furthermore, since the current LPP model does not include either droplet
formation or droplet-interface interaction, only droplets that are at a given
distance away from the liquid jet interface will be converted. The distance is
typically chosen to be the droplet diameter $d_{p}$. The overall conversion
criteria can then be expressed as

$\displaystyle\\{\mathrm{LPP\
Qualified}\\}=\\{V_{p}<V_{crit}\\}\,\,\&\&\,\,\\{|\gamma_{p}-1|<\epsilon_{\gamma}\\}\,\,\&\&\,\,\\{{\bf
x}_{p}\in\mathcal{R}_{ai}\\}\,,$ (74)

where $\epsilon_{\gamma}$ is the tolerance for the aspect ratio, and
$\mathcal{R}_{ai}$ is the region away from the interface.

##### RD-to-LPP conversion

If a RD satisfies the LPP criteria, the volume fraction $C$ in the
corresponding cells is set to zero. A new LPP droplet is created and added to
the LPP array, together with its volume and velocity. Furthermore, the flow
field in the same region will be replaced by the undisturbed flow, so that the
LPP droplet sees the proper undisturbed velocity field. This is done by
interpolating each component of the fluid velocity along the corresponding
coordinate from the surface into the interior of a cubic region centered at
the droplet location. The reconstructed velocity field is globally divergence-
free if the velocity on the surface of the interpolation region is divergence-
free. Unless the particle Reynolds number is very small, a cubic region with
an edge size twice the droplet diameter is sufficient to reconstruct the
undisturbed flow field in a satisfactory way.

##### LPP-to-RD conversion

When a LPP droplet does not satisfy the LPP criteria anymore, for example when
it is too close to the liquid-jet interface, it will be converted back to a
resolved droplet. First, a spherical droplet is rebuilt around ${\bf x}_{p}$,
by specifying the volume fraction $C$ in the cells that will be occupied by
the resolved droplet. The velocity field needs to be updated in the same
cells. To account for the momentum of both the droplet and the virtual fluid
moving with the droplet, the velocity in the cells occupied by the resolved
droplet should be changed to ${\bf u}^{\prime}_{p}$, which is calculated as

${\bf u}^{\prime}_{p}-{\bf
u}=\Big{(}1+\alpha\,\frac{\rho}{\rho_{p}}\Big{)}\,({\bf u}_{p}-{\bf u})\,,$
(75)

where $\alpha$ is the ratio between the volumes of the virtual fluid and the
droplet. For the added-mass effect, it is considered that the ratio between
the volumes of the virtual fluid to be accelerated through the inviscid
mechanism and the particle is the added-mass coefficient $C_{M}$, which is
equal to 0.5 for a sphere in incompressible flows. Due to finite viscous
diffusion time scale, the history force usually needs to be expressed in
integral form. Nevertheless, it has been shown in [49] that, if the particle
and ambient fluid acceleration time scales are much larger than the viscous
diffusion time scale, the history force can also be expressed as a non-
integral form like the added-mass force, and a viscous-unsteady coefficient
$C_{vu}$ similar to the added-mass coefficient $C_{M}$ can be derived as

$\displaystyle C_{vu}\approx
8.51\left(\frac{0.75+0.105{\mathrm{Re}_{p}}}{{\mathrm{Re}_{p}}}\right)\,.$
(76)

The viscous-unsteady coefficient $C_{vu}$, can be considered as the ratio
between the volumes of the virtual fluid to be accelerated through the viscous
mechanism and the particle. The excess momentum added through ${\bf
u}^{\prime}_{p}$ is to mimic the effects of the added-mass and history forces
on accelerating the surrounding fluid around the droplet. Therefore, $\alpha$
can be estimated as the sum of $C_{M}$ and $C_{vu}$ and thus depends on the
droplet Reynolds number as well.

##  4 Testing

The testing of the code is performed automatically. The short version of
testing is done by typing make test, the long version by typing make longtest.
The short version takes approximately 5 minutes on a laptop with an i7
processor and the long version takes approximately 25 minutes. All the
resulting tests give a report of “PASS” or “FAIL”. Each test is contained in a
subdirectory of the Test directory. The subdirectory corresponding to a test
has a self-explanatory name, e.g. PresPoiseuille for the pressure-driven
Poiseuille flow.

The tests can be divided into two categories: elementary tests which are
basically sanity checks verifying that the code is not corrupted and finds
elementary flows easily, and more complex test that are in some cases a
verification of the code, comparing it to analytical solutions. However, the
verification has not been pushed very far, since the code is an assembly of
methods that have been tested extensively elsewhere, see for example TSZ for a
review of these tests. The more recent methods such as the “mass-momentum
consistent” option for velocity advection, has been tested extensively in [9]
although several test cases of the method are included in the test suite and
will be described below. In order to avoid the introduction of more equations
we assume a constant surface tension coefficient $\sigma$ in all the testing
section.

In many of the tests the solution computed during the test run is compared to
a reference solution. In some cases the reference solution was computed by the
authors at a different resolution and is included with the code distribution.
In other cases the “reference” is an near-identical numerical solution
performed by the authors and included in the code distribution. In that case
the reference should be identical to the solution except in some cases where
tiny changes in the code or the implementation create moderately large
differences. This is the case for the Raindrop test below.

###  4.1 Elementary tests

####  4.1.1 Poiseuille flow

An elementary Poiseuille flow [50] is tested. The simulation is set up in a
$8\times 8\times 2$ grid in which the system reduces to a 2D planar flow in
the box $(0,1)\times(0,1)$. The parameters are $\|\nabla p\|=\mu=\rho=1$. The
boundary conditions set pressure on the left and right. The simulation is
continued until the flow becomes stationary, which happens with $10^{-3}$
accuracy around time $1$. This time is reached in $1000$ time steps. The
explicit version of viscous diffusion is used.

This test passes if the numerical segments are tangent to the theoretical
profile as seen in Fig. 13. Because the Poiseuille flow profile is quadratic,
second-order finite differences offer exact values for the second derivative
of velocity, which ensures that the profile found is exactly a parabola. The
accuracy of the parabola amplitude is set by the quality of the approximation
of the $u=0$ boundary condition (solid wall). Since this condition is set at
first order, there is a small ${\cal O}(h^{2})$ difference with the exact
parabola.

The bottom wall tangential velocity boundary condition is $u=0$ on $y=0$.
Since the wall is at $y_{i,1/2}$, the boundary condition is written using a
ghost velocity at $y_{i,-1/2}$ that satisfies $u^{\rm
ghost}_{i,-1/2}+u_{i,1/2}=2u^{\rm wall}$. The boundary condition is thus
implemented by writing in a “ghost layer” of the grid

$u^{\rm ghost}_{i,-1/2}\leftarrow 2u^{\rm wall}-u_{i,1/2}.$ (77)

The result is shown on Fig. 13.

Figure 13: Poiseuille flow test case.

It is also possible to run the Poiseuille flow test in other manners, for
example with set inflow velocity, or to run it with the implicit version of
viscous diffusion, in which case the flow converges in a few time steps.

####  4.1.2 Stokes flow around a disk

A pressure driven flow around a disk of diameter $1/2$, with the other
parameters as before, is set and left to evolve until steady state. The
advection operator ${\cal L}_{\rm adv}$ is turned off which ensures that the
steady state is a Stokes flow. The explicit version of the viscous terms is
used. The test resides in the test directory PresDisk. If the implicit version
is used, convergence to the steady state can be faster but an error of order
$\tau$ affects the solution. The result is shown on Fig. 14.

Figure 14: Stokes flow around a disk test case.

There are two other similar tests beyond this first one. In the second one,
located in the directory Disk, the flow is driven by a body force akin to
gravity instead of being driven by pressure. This second test is still with
periodic boundary conditions. The third test has inflow and outflow conditions
and is located in the test directory Inflowdisk.

####  4.1.3 Droplet advection

This is an elementary test to check whether the VOF method is operating
normally. The final state of the volume fraction field $C$ is compared to a
precomputed value. The test has to be visualized “by hand” by the user, with
the help of graphics software such as Visit or Paraview. One should see an
undeformed droplet moving across the domain.

####  4.1.4 Cylinder advection

This is a more sophisticated test that probes the “mass-momentum consistent”
option. The test is described in detail in [9]. If the velocity field stays
uniform and constant and as a result the droplet is advected undeformed, it
means that momentum $\rho{\bf u}$ and density $\rho$ are advected in lock-
step, so that the operation $\rho{\bf u}/\rho$, at each time step, gives the
constant ${\bf u}$.

####  4.1.5 Speed

This is not so much a test as a report on the code speed. On a 2015 MacBook, a
single processor run delivered a speed of $1.6\,10^{6}$ cells per second. On
an AMD EPYC 7351, a single processor run delivered a speed of $1.9\,10^{6}$
cells per second and a parallel eight-processor run delivers a speed of
$1.65\,10^{6}$ cells per second and per processor (The parallel test is
optional and can be run by the user by typing sh run-speed-test.sh $N$ where
$np=N^{3}$ is the number of processors desired). A significant drop of the
code speed from this value would be a serious issue.

###  4.2 Capillary wave

In this section we present results of the oscillation of planar capillary
waves between two viscous fluids with equal density and viscosity in the
presence of surface tension. The interface between the two fluids is slightly
perturbed with a sinusoidal function of small amplitude $a_{0}$ and the
initial velocity is set to zero. The solution of this problem is governed by
the Laplace number which is $La=\sigma\rho\lambda/\mu^{2}=3000$, where
$\lambda$ is the wavelength of the sinusoidal function. Simulation are
performed in a box of dimensions $L_{x}=\lambda$ and $L_{y}$. The results are
compared to the analytical initial-value solution (AIVS) obtained in [51, 52]
for small-amplitude capillary waves in viscous fluids. In the AIVS the aspect
ratio $L_{y}/L_{x}$ should be sufficiently large (at least a value of 2) and
the initial capillary wave amplitude $a_{0}$ sufficiently small. Moreover the
time step $\tau$ and the tolerance of the solvers should be at convergence. We
have checked that all these four parameters were at convergence for fixed $h$.
We then look at the dependence of the remaining error on $h$. Figure 15
compares the temporal evolution of the amplitude of the interface perturbation
with the AIVS.

Figure 15: Comparison of the temporal evolution of the normalized amplitude of the capillary wave obtained numerically and of the analytical initial-value solution by [51]. $N=\lambda/h$ | 8 | 16 | 32 | 64 | 128
---|---|---|---|---|---
Gerris | 0.159 | 0.032 | 0.0077 | 0.0022 | $5.5\,10^{-4}$
Basilisk | 0.139 | 0.024 | 0.0069 | 0.0024 | $4.8\,10^{-4}$
Paris | 0.050 | 0.023 | 0.0054 | 0.0015 | $4.1\,10^{-4}$
Table 1: Relative $L_{2}$ error of the numerical solution for capillary waves
for various codes and numbers of grid points $N$ per wavelength. The errors
estimated by the codes have been rounded to the nearest digit. Results for
Gerris are from the website http://gfs.sf.net, not from the paper [34].
Results for Basilisk, have been obtained by the authors from the code
published on the website http://basilisk.fr.  Figure 16: Relative $L_{2}$
error of the numerical solution as a function of the number of grid points per
wavelength $N=\lambda/h$ for the capillary wave test.

The relative $L_{2}$ error norm of the difference between the numerical
solution of a few numerical codes and the AIVS has also been computed. The
results, depicted in Table 1 and Fig. 16, show second-order convergence for
coarse grids. For very refined grids, when the error is below 1%, the accuracy
on the solution is controlled by the initial amplitude of the wave. In the
case of the resolution $\lambda/h=128$, instead of a $0.01$ amplitude as in
[34], a 0.005 amplitude had to be used to match the AIVS.

###  4.3 Oscillating droplets and bubbles

$D$ | $\mu_{l}$ | $\mu_{g}$ | $\rho_{l}$ | $\rho_{g}$ | $\sigma$
---|---|---|---|---|---
$(m)$ | $(kg\,m^{-1}s^{-1})$ | $(kg\,m^{-1}s^{-1})$ | $(kg\,m^{-3})$ | $(kg\,m^{-3})$ | $(kg\,s^{-2})$
$3\,\,10^{-3}$ | $10^{-3}$ | $1.7\,\,10^{-5}$ | $10^{3}$ | $1.2$ | 0.0728
Table 2: Physical parameters (defined in the text) for the oscillating droplet.  $r$ | $m$ | ${\rm La}\,$
---|---|---
$\rho_{l}/\rho_{g}$ | $\mu_{l}/\mu_{g}$ | $\sigma\rho_{l}d/\mu_{l}^{2}$
833.3 | 58.82 | 218400
Table 3: Dimensionless parameters for the oscillating droplet, ${\rm La}\,$ is
the Laplace number.

A droplet with a large density ratio is initialized with a small ellipsoidal
deformation. The droplet has air-water properties described in Tables 2 and 3.
The initial shape, when tracked with the Front, is shown in Fig. 17. This test
is not designed to assess the accuracy of the code, since the implemented
numerical methods have already been used and tested elsewhere (see for example
[53] for an oscillating droplet test with the VOF method and [54, 55] for
similar tests with Front Tracking. However, it should be noted that in these
references the tests are 2D, hence easier). We thus expect the accuracy to be
similar to that of already published and tested codes using similar curvature
and surface tension methods. The purpose of the test is rather to ensure that
the code is working as expected, and to assess whether air-water properties,
which are often creating stability problems, are in fact in the stable regime
of the code.

Figure 17: Initial ellipsoidal shape of the droplet or bubble with the Front.
Figure 18: Amplitude of capillary oscillations of the Droplet test case.
Figure 19: Amplitude of capillary oscillations of the FrontDroplet test case.

Figure 18 shows the amplitude oscillations as a function of time for a droplet
of $D/h=19.2$ grid points per diameter and an initial excentricity of $0.75$.
The test is run without the mass-momentum consistent option, which is here
seen to be unnecessary for the stability, and with the mixed-height option
discussed in A.3. The reference solution, plotted alongside the test
simulation result, is obtained from the same VOF simulation at the larger
resolution $D/h=38.4$. Satisfying agreement is found.

Note that when this test is run automatically in the test suite, the reference
solution stored in the Test/Droplet directory has been obtained by our code
running in the same conditions, a device frequently used in several test cases
in the code test suite. That way, one tests that the behavior of the code has
not changed drastically, but not that the code (original version and current
version) is correct.

We test the Front-Tracking part of the code by simulating the same droplet
test. Results are shown on Fig. 19. The reference solution is once again the
VOF simulation at larger resolution ($D/h=38.4$). Satisfying agreement is
found in this test as well.

We repeat these tests by just inverting the phases, thus initializing an air
bubble inside water. The physical and the numerical parameters (such as the
scheme options) are the same as before.

Figure 20: Kinetic energy of capillary oscillations of the Bubble test case.
Figure 21: Kinetic energy of capillary oscillations of the FrontBubble test
case.

Figure 20 shows for the Bubble test case the kinetic energy oscillations as a
function of time for a bubble of $D/h=19.2$ grid points per diameter and an
initial excentricity of $0.75$. The kinetic energy is used this time instead
of the deformation amplitude. The reference solution is again obtained from
the VOF simulation at larger resolution, with $D/h=38.4$ grid points. There is
again good agreement.

We test again the Front-Tracking part of the code by simulating the same
bubble test. Results are shown on Fig. 21. The reference solution is again
that of the previous case (VOF simulation with resolution $D/h=38.4$). The
agreement is satisfying, but a small drift of the kinetic energy is observed.

###  4.4 Falling raindrop

This is an important test case, since it is very demanding for several codes.
A spherical raindrop is setup with a diameter $d=3\,$mm and allowed to fall in
air under gravity. The droplet should remain approximately spherical, with at
most a pancake or bun-like shape, but for many numerical codes this case is
rather difficult and spurious atomisation of the droplet is seen. For an
example with the Basilisk code see [56]. For details about the setup of the
test we refer the reader to [9]. The test is performed with the parameters of
Table 4. Notice that these parameters differ slightly from those of [9] as we
input here the values of the physical parameters for air and water at
temperature $T=20\,^{\circ}C$. The grid resolution is low $D/\Delta x=8$,
since the test is more demanding (i.e., it leads to a blowup of the code more
easily) at low resolutions. Figure 22 provides the result of the Raindrop test
in the code distribution. It is seen that the solution deviates somewhat from
the reference, however this is not worrisome, as explained in the introduction
to this section 4, since the flow is in a regime that is very sensitive to
physical parameters and initial conditions and any small change in the code
will create a deviation of that sort. The shape of the droplet is shown in
Fig. 23.

$\mu_{l}$ | $\mu_{g}$ | $\rho_{l}$ | $\rho_{g}$
---|---|---|---
$(kg\,m^{-1}s^{-1})$ | $(kg\,m^{-1}s^{-1})$ | $(kg\,m^{-3})$ | $(kg\,m^{-3})$
$1.0016\,10^{-3}$ | $1.835\,\,10^{-5}$ | $998.2$ | $1.19$
Table 4: Physical parameters for the raindrop test. Only the parameters that
differ from Table 2 are given.  Figure 22: Kinetic energy of a $3\,$mm falling
raindrop at low resolution $D/\Delta x=8$. Figure 23: Simulated shape of a
$3\,$mm falling raindrop at low resolution $D/\Delta x=8$.

##  5 Installation and Usage

###  5.1 External Libraries

When solving elliptic equations, we may apply the Hypre package SMG [57], a
semi-coarsening multigrid solver with 3D plane smoothing on structured, cuboid
meshes, as mentioned in Section 3.7.2.

Installation of the static library versions (*.a files) of Hypre is controlled
via the Makefile. For example, for a Linux system, this can be realized via
the line

    HYPRE_DIR = $(HOME)/some_path/hypre-2.10.0b/src/lib

in the Makefile. Note that the file libHYPRE.a is placed in the directory
HYPRE_DIR. Consequently, the actual linking is ensured by another Makefile
line

    HYPRE_LIBS = -L$(HYPRE_DIR) -lHYPRE

which is specified in the default Paris distribution. Note that the Hypre
version used for the majority of Paris development has been 2.10.1; it has
been tested for stability in serial and massively parallel runs [42]. In case
of Hypre-related problems, fallback to version 2.10 is recommended for
debugging.

The Vofi library [29, 30] may be used to initialize the volume fraction field
$C$ in 2D and in 3D (see Sec. 3.5.3). It is an interesting option in cases
where initial conditions depend strongly on a very precise interface geometry,
e.g. free-surface solutions of bubble dynamics [43, 42], and more rarely when
initialization requires a considerable amount of computational time compared
to Paris simulation time (as in the curvature test case). Linking of this
library is performed in the exact same fashion as in the case of Hypre. The
static library file is named libvofi.a, and the respective linker switch is
-lvofi.

Since Hypre and Vofi are indeed optional, compilation without them is made
easy. The user may set or unset the variables HAVE_VOFI and HAVE_HYPRE in his
shell prior to compiling. If these variables are set the Makefile will attempt
to locate the corresponding libraries, if they are not set the compilation
proceeds without the libraries. The fallback for Hypre is the built-in Gauss-
Seidel solver followed by the in-code multigrid solver, and the fallback for
Vofi are the built-in VOF initialization procedures.

###  5.2 Output Files

Various output formats are available in the code: VTK, SILO, and MPI I/O.
While the VTK option generates ASCII files, the latter two produce binary
files. The output in SILO format is done based on the SILO library developed
at LLNL. For both the VTK and SILO output options, the independent parallel
approach is used, namely every MPI process generates a separate file. This
will become an issue for large-scale simulations using many MPI processes,
since creating a large number of small files simultaneously may crash the
indexing server. An output option based on the MPI I/O library is implemented
in PARIS for large-scale simulations, which adopts, instead, a cooperative
parallel approach and creates a single file for each variable for each output.
A post-processing code was developed in PARIS to convert the MPI I/O outputs
and SILO files offline for visualization. This code is available in the
distribution as the file util/post_utility.f90.

###  5.3 Input Files

Paris requires a set of input files (in text format) to initialize and start
the simulation. These files are

  * 1.

input \- global and front-solver parameters.

  * 2.

inputFS \- free-surface solver parameters (optional, free-surface simulations
only).

  * 3.

inputlpp \- Lagrangian point-particle module parameters (optional, implicitly
requires two-phase flow, see [44]).

  * 4.

inputvof (optional, VOF parameters, as above).

  * 5.

inputsolids (optional, solid objects parameters).

All input files contain over 220 parameters, thus listing all of them is not
practical in this paper; instead we will only point to general rules governing
the use of these files. The reader may find the default values of these
parameter in the source code, usually just below the namelist instruction.

All the input lines contain the parameter specifications written as “variable
= value”, with values being reals, integers, string or boolean (T/F). In most
cases (in the code versions distributed in main darcs tree) the variables are
commented Fortran-style, i.e. in the same line, following an exclamation mark,
for example

    MaxDt   = 5.e-5  ! maximum size of time step
    dtFlag  = 2      ! 1: fixed dt; 2: fixed CFL
    dt      = 1.0e-4 ! dt in case of dtFlag=1
    MAXERROR= 1.0d-6 ! Residual for Poisson solver
    CFL     = 0.042Ψ

It must be noted that the Paris source code uses the Fortran namelist input
type, consequences of that being that all input files have “free format”, i.e.
lines can change order or be deleted. All variables are initialized to default
values in the source code. Thus, Paris will initialize with an empty input
file, however in such a case the simulation will be short. Indeed by default
Nx=0 (the grid has zero points in $x$ direction) in order to prevent
simulations with some of the absurd input files that could be selected by
mistake. Thus a minimum input file should contain at least a specification of
Nx. For more demanding simulations, beginner users are encouraged to
familiarize themselves with input file examples, such as templates found in
the Tests sub-directory in the distribution which can be copied and modified
to create new Paris cases. Note that in the Test suite, input files are often
generated from template files such as inputfilename.template using shell
scripts.

##  6 Acknowledgements

We thank Dr. V. Le Chenadec, Mr. C. Pairetti, Dr. S. Popinet and Dr. S.
Vincent for useful conversations on the topics of this paper.

Portions of this work were supported by National Science Foundation Grants
CBET-1335913 and NSF DMS-1620158, by the ANR MODEMI project (ANR-11-MONU-0011)
program and by 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 like to
acknowledge the MESU computing facilities of Sorbonne Université. Finally, the
simulation data are visualized by the software VisIt developed at the Lawrence
Livermore National Laboratory.

## References

  * [1] R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Annu. Rev. Fluid Mech. 31 (1999) 567–603.
  * [2] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct Numerical Simulations of Gas-Liquid Multiphase Flows, Cambridge University Press, 2011.
  * [3] S. Dabiri, D. Fuster, Y. S. Ling, L. Malan, R. Scardovelli, G. Tryggvason, P. Yecko, S. Zaleski, PARIS Simulator Code: A PArallel Robust Interface Simulator that combines VOF and Front-Tracking, http://parissimulator.sf.net (2012-2015).
  * [4] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, G. Zanetti, Modelling merging and fragmentation in multiphase flows with SURFER, J. Comput. Phys. 113 (1994) 134–147.
  * [5] S. Popinet, The gerris flow solver, http://gfs.sf.net (2001-2014).
URL http://gfs.sf.net

  * [6] S. Popinet, Basilisk, a Free-Software program for the solution of partial differential equations on adaptive Cartesian meshes (2018).
URL http://basilisk.fr

  * [7] J. Brackbill, D. B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335–354.
  * [8] R. P. Fedkiw, T. Aslam, B. Merriman, S. Osher, A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method), Journal of computational physics 152 (2) (1999) 457–492.
  * [9] T. Arrufat, M. Crialesi-Esposito, D. Fuster, Y. Ling, L. Malan, S. Pal, R. Scardovelli, G. Tryggvason, S. Zaleski, A mass-momentum consistent, volume-of-fluid method for incompressible flow on staggered grids, Computers & Fluids 215 (2021) 104785.
  * [10] 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.
  * [11] O. Desjardins, V. Moureau, Methods for multiphase flows with high density ratio, Center for Turbulent Research, Summer Programm 2010 (2010) 313–322.
  * [12] 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.
  * [13] 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.
  * [14] S. Ghods, M. Herrmann, A consistent rescaled momentum transport method for simulating large density ratio incompressible multiphase flows using level set methods, Physica Scripta 2013 (T155) (2013) 014050.
  * [15] 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.
  * [16] J. K. Patel, G. Natarajan, A novel consistent and well-balanced algorithm for simulations of multiphase flows on unstructured grids, Journal of computational physics 350 (2017) 207–236.
  * [17] N. Nangia, B. E. Griffith, N. A. Patankar, A. P. S. Bhalla, A robust incompressible navier-stokes solver for high density ratio multiphase flows, Journal of Computational Physics.
  * [18] C. B. Ivey, P. Moin, Conservative and bounded volume-of-fluid advection on unstructured grids, Journal of Computational Physics 350 (2017) 387–419.  doi:10.1016/j.jcp.2017.08.054.
  * [19] M. Owkes, O. Desjardins, A mass and momentum conserving unsplit semi-Lagrangian framework for simulating multiphase flows, Journal of Computational Physics 332 (2017) 21–46.  doi:10.1016/j.jcp.2016.11.046.
  * [20] M. Rudman, A volume-tracking method for incompressible multi-fluid flows with large density variations, Int. J. Numer. Meth. Fluids 28 (1998) 357–378.
  * [21] 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.
  * [22] M. Razizadeh, S. Mortazavi, H. Shahin, Drop breakup and drop pair coalescence using front-tracking method in three dimensions, Acta Mechanica 229 (2018) 1021–1043.
  * [23] J. Lu, G. Tryggvason, Direct numerical simulations of multifluid flows in a vertical channel undergoing topology changes, Physical Review Fluids 3 (2018) 084401 (20 pages).
  * [24] J. Lu, G. Tryggvason, Multifluid flows in a vertical channel undergoing topology changes–effect of void fraction, Physical Review Fluids 4 (2019) 084301\.  doi:https://doi.org/10.1103/PhysRevFluids.4.084301.
  * [25] D. M. McQueen, C. S. Peskin, A three-dimensional computational method for blood flow in the heart: (ii) contractile fibers, Journal of Computational Physics 82 (1989) 289–297.
  * [26] E. Aulisa, S. Manservisi, R. Scardovelli, S. Zaleski, Interface reconstruction with least-squares fit and split advection in three-dimensional cartesian geometry, J. Comput. Phys. 225 (2007) 2301–2319.
  * [27] D. L. Youngs, An interface tracking method for a 3D Eulerian hydrodynamics code, Tech. rep., AWRE, technical Report 44/92/35 (1984).
  * [28] R. Scardovelli, S. Zaleski, Analytical relations connecting linear interfaces and volume fractions in rectangular grids, J. Comput. Phys. 164 (2000) 228–237.
  * [29] 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.
  * [30] 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.
  * [31] 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.
  * [32] R. Scardovelli, S. Zaleski, Interface reconstruction with least-square fit and split Lagrangian-Eulerian advection, Int J. Numer. Meth. Fluids 41 (2003) 251–274.
  * [33] S. Popinet, Numerical models of surface tension, Annual Review of Fluid Mechanics 50 (2018) 49–75.
  * [34] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, J. Comput. Phys. 228 (2009) 5838–5866.
  * [35] G. Bornia, A. Cervone, S. Manservisi, R. Scardovelli, S. Zaleski, On the properties and limitations of the height function method in two-dimensional cartesian geometry, J. Comput. Phys. 230 (2011) 851–862.
  * [36] M. Owkes, O. Desjardins, A mesh-decoupled height function method for computing interface curvature, J. Comput. Phys. 281 (2015) 285 – 300.  doi:http://dx.doi.org/10.1016/j.jcp.2014.10.036.
URL http://www.sciencedirect.com/science/article/pii/S0021999114007189

  * [37] W. L. Briggs, A multigrid tutorial, SIAM Philadelphia, 1987.
  * [38] M. Kang, R. P. Fedkiw, X.-D. Liu, A boundary condition capturing method for multiphase incompressible flow, Journal of Scientific Computing 15 (3) (2000) 323–360.
  * [39] R. K. C. Chan, R. L. Street, A computer study of finite-amplitude water waves, J. Comput. Phys. 6 (1970) 68–94.
  * [40] S. Popinet, S. Zaleski, Bubble collapse near a solid boundary: a numerical study of the influence of viscosity, J. Fluid Mech. 464 (2002) 137–163.
  * [41] M. Sussman, A second order coupled level set and volume-of-fluid method for computing growth and collapse of vapor bubbles, J. Comput. Phys. 187 (2003) 110–136.
  * [42] W. Aniszewski, S. Zaleski, A. Llor, L. Malan, Numerical simulations of pore isolation and competition in idealized micro-spall process, International Journal of Multiphase Flowdoi:https://doi.org/10.1016/j.ijmultiphaseflow.2018.10.013.
URL http://www.sciencedirect.com/science/article/pii/S0301932218303082

  * [43] L. Malan, Y. Ling, R. Scardovelli, A. Llor, S. Zaleski, Detailed numerical simulations of pore competition in idealized micro-spall using the VOF method, Computers & Fluids 189 (2019) 60–72.
URL arXiv:1711.04561[physics.flu-dyn]

  * [44] 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

  * [45] M. Herrmann, A parallel Eulerian interface tracking/Lagrangian point particle multi-scale coupling procedure, J. Comput. Phys. 229 (3) (2010) 745–759.
  * [46] C. T. Crowe, M. Sommerfield, Y. Tsuji, Multiphase flows with droplets and particles, CRC Press, 1998.
  * [47] R. Clift, W. H. Gauvin, The motion of particles in turbulent gas streams, Proc. Chemeca 1 (1970) 14–28.
  * [48] M. R. Maxey, B. K. Patel, E. J. Chang, L. P. Wang, Simulations of dispersed turbulent multiphase flow, Fluid Dyn. Res. 20 (1997) 143–156.
  * [49] Y. Ling, M. Parmar, S. Balachandar, A scaling analysis of added-mass and history forces and their coupling in dispersed multiphase flows, Int. J. Multiphase Flow 57 (2013) 102–114.
  * [50] P. Kundu, I. Cohen, D. Dowling, Fluid Mechanics, 891 pp, Elsevier, Amsterdam, 2012\.
  * [51] A. Prosperetti, Motion of two superposed viscous fluids, Phys. Fluids 24 (1981) 1217–1223.
  * [52] F. Denner, G. Paré, S. Zaleski, Dispersion and viscous attenuation of capillary waves with finite amplitude, The European Physical Journal Special Topics 226 (6) (2017) 1229–1238.
  * [53] D. Fuster, G. Agbaglah, C. Josserand, S. Popinet, S. Zaleski, Numerical simulation of droplets, bubbles and waves: state of the art, Fluid Dyn. Res. 41 (6) (2009) 065001.
  * [54] D. J. Torres, J. U. Brackbill, The Point-Set Method: Front-Tracking without Connectivity, Journal of Computational Physics 165 (2) (2000) 620–644.
  * [55] U. Olgac, D. Izbassarov, M. Muradoglu, Direct numerical simulation of an oscillating droplet in partial contact with a substrate, Computers and Fluids 77 (C) (2013) 152–158.
  * [56] C. Pairetti, S. Popinet, S. M. Damián, N. Nigro, S. Zaleski, Bag mode breakup simulations of a single liquid droplet, 6th European Conference on Computational Mechanics (ECCM 6) and 7th European Conference on Computational Fluid Dynamics (ECFD 7) 1115 June 2018, Glasgow, UK.
  * [57] P. Sloot, C. Tan, J. Dongara, A. H. (Editors), Hypre: a library of high performance preconditioners, in: Computational Science - ICCS, Springer-Verlag, ICSS 2002, Berlin, 2002.

##  Appendix A Details of curvature computation from height functions (HF)

###  A.1 Height computation

The height computation is performed as described in Section 3.6.2. A more
general definition than (50) is to consider for each cell $\Omega_{i,j,k}$ the
possible existence of one of six height functions defined with reference to a
direction ${\bf e}_{a}$, $1\leq a\leq 3$, where ${\bf e}_{a}$ is one of the
Cartesian base vectors aligned with the grid, and an orientation
$\epsilon=-{\rm sign}({\bf e}_{a}\cdot\nabla C)$ (the minus sign ensures that
the canonical situation where the “liquid” $C=1$ is below the “air” $C=0$ has
$\epsilon=1$. It also corresponds to the sign convention for the interface
normal ${\bf n}\delta_{S}=-\nabla\chi$). This cell-and-orientation-dependent
HF is defined as

$H_{i,j,k}^{(a,+)}=\sum_{{\rm stencil}\,S}C_{l,m,n}-L_{i,j,k}$ (78)

where the sum is over all the cells in a one-dimensional symmetric “stencil”
or “stack” $S$ centered on ${\bf x}_{i,j,k}$ and oriented parallel to ${\bf
e}_{a}$ and for the “positive” orientation $\epsilon$.
$L_{i,j,k}=\epsilon({\bf x}_{i,j,k}-{\bf x}_{O})\cdot{\bf e}_{a}$ is the
distance from the base of the stack to the cell center ${\bf x}_{i,j,k}$. When
the orientation is $\epsilon=-1$

$H_{i,j,k}^{(a,-)}=\sum_{{\rm stencil}\,S}\left(1-C_{l,m,n}\right)-L_{i,j,k}$
(79)

and the distance $L_{i,j,k}$ is now with reference to an origin in direction
$-{\bf e}_{a}$ from the cell. An example of stack is shown on Fig. 6(c). This
height function can be computed whenever both the bottom cell and the top cell
of the stack do not contain the interface, and the interface crosses only once
the intermediate cells. This can be tested by requiring that there are cells
with $C=0$ and $C=1$ at the top and bottom of the stack (see again Fig. 6(c)).

For a straight line interface the height function is exact. It is interesting
to see how many cells are needed in the stack $S$ to find the height for a
straight line in 2D. The most “difficult” case is the $45^{\circ}$ case. Thus,
considering a cell crossed by the interface, one should explore one cell above
and one cell below that cell. With the addition of the full and empty cells
this requires the exploration of two cells above and below the starting cell.
The total number of cells for a symmetric stencil about the starting cell
would thus be five, but the total number of cells in a stencil maybe as low as
three. However, even with a vanishing amount of curvature five cells in a
symmetric stencil are not sufficient and seven cells are needed. Thus one
would need to explore $N_{d}=3$ cells above and below. In three dimensions the
most “dangerous” cases now have the plane normal ${\bf n}=(1,1,1)$. A similar
reasoning also leads to a height of seven cells for the symmetric stencil in
3D. Note that the stencil does not have to be symmetric, rather this is a
accidental feature or our implemention of the method.

For each cell, it is determined whether there is a full or empty cell at a
maximum distance $N_{d}$ (the parameter NDEPTH in the code) above or below the
cell. The value NDEPTH=3 is hard-coded. In order to function in parallel, and
since only two layers of cells are exchanged between MPI processes, the sum in
(50) is broken in two parts, one in each processor. Then the processes
exchange two pieces of information, the “partial sums” just computed and the
lengths $L_{i,j,k}$ in expression (78), allowing them to reconstruct the full
sum.

###  A.2 First pass, first attempt: fully-aligned heights

In the first pass, a loop is performed over all cells $\Omega_{i,j,k}$ cut by
the interface, hence having $0<C_{i,j,k}<1$. In the first pass two attempts
are made. Let the current cell be $\Omega_{0}$ with grid coordinates
$i_{0},j_{0},k_{0}$. In the first attempt, the normal ${\bf n}$ is estimated
by MYCS [26]. Then the grid direction ${\bf e}_{a}$ closest to the normal is
determined (maximum of $|{\bf n}\cdot{\bf e}_{a}|$). Without loss of
generality we consider the case $a=3$ and the horizontal plane perpendicular
to ${\bf e}_{3}$, that is the grid plane most closely aligned with the
interface. A $3\times 3$ planar block of cells $\Sigma_{0}$, aligned with this
plane, is selected containing the nine cells $\Omega_{i,j,k_{0}}$, with
$i_{0}-1\leq i\leq i_{0}+1$ and $j_{0}-1\leq j\leq j_{0}+1$. For all these
cells, either a height $H_{i,j}$ is readily available, or is searched in the
above and below cells over two layers, that is for $k_{0}\pm 1$ and $k_{0}\pm
2$. When all nine heights are available, the following coefficients of the
polynomial (51) can be found using finite differences

$\displaystyle a_{1}$ $\displaystyle=\partial^{2}_{xx}H\simeq
H_{1,0}-2H_{0,0}+H_{-1,0}\,,$ $\displaystyle a_{2}$
$\displaystyle=\partial^{2}_{yy}H\simeq H_{0,1}-2H_{0,0}+H_{0,-1}\,,$
$\displaystyle a_{3}$
$\displaystyle=\partial^{2}_{xy}H\simeq(H_{1,1}-H_{-1,1}-H_{1,-1}+H_{-1,-1})/4\,,$
$\displaystyle a_{4}$
$\displaystyle=\partial_{x}H\simeq(H_{1,0}-H_{-1,0})/2\,,$ $\displaystyle
a_{5}$ $\displaystyle=\partial_{y}H\simeq(H_{0,1}-H_{0,-1})/2\,.$ (80)

###  A.3 First pass, second attempt: mixed heights

If the first attempt fails then a “mixed heights” approach is used, but only
if the parameter MIXED_HEIGHTS is set to ’T’. For every height
$H_{i,j,k}^{(a,\epsilon)}$ that has been calculated, a point ${\bf
x}_{i,j,k}^{a}=H_{i,j,k}^{(a,\epsilon)}{\bf e}_{a}+x_{b}{\bf e}_{b}+x_{c}{\bf
e}_{c}$ is defined, where $x_{b}$ and $x_{c}$ are the cell-central coordinates
in the two directions other than ${\bf e}_{a}$. Since there are six height
directions, up to 54 points could be computed. However, a general orientation
is computed using the MYCS normal and only the orientations compatible with
that orientation are retained, which yields 27 possible points. With certain
point configurations there is a risk of a degenerate case where the least-
square linear operator is not invertible. This happens for example in the set
of six points obtained with combinations of $x=0,1$, $y=-1,0,1$, $z=0$. All
paraboloids of the form $z=x(1-x)$ pass through these points. Other
degeneracies are possible: points all on a circle will be fitted by infinitely
many revolution paraboloids $z^{\prime}=\kappa(x^{2}+y^{2})/2$. To avoid these
degeneracies, and after trial and error, the minimum number of points
requested is hard-coded as $N_{s}=$NFOUND_MIN +1$=7$. In addition to these
“mixed heights” points, the centroid of the VOF face in the current cell
$\Omega_{i,j,k_{0}}$ is added to the set of points to be fitted. In some
cases, different directions ${\bf e}_{a}$ could yield two close points, say
${\bf x}_{i,j,k}^{a}$ and ${\bf x}_{l,m,n}^{b}$, in the same cell or in a
neighboring cell. In this case, if $||{\bf x}_{i,j,k}^{a}-{\bf
x}_{l,m,n}^{b}||<h/2$ then one of the two points is rejected. Which point is
rejected depends on the order in which points are added to the stack, which in
turn depends on the order in which mixed heights are investigated, typically
closest to the general orientation. Before the fitting is performed, an
approximate normal is computed using the MYCS approach and the coordinate
system is rotated so that the $z$ axis is now aligned with the approximate
normal. The rotation is optional and is controlled by the parameter
DO_ROTATION. We have found that performing rotation has a certain positive
influence on the accuracy of the results, although it is not clear why.

By default the parameter MIXED_HEIGHTS is set to ’T’ (true). This gives less
accurate results in $L_{1}$ norm for curvature, but a smoother computation and
as a result simulations appear to be more stable for large density ratios when
the mass-momentum consistent scheme is not used. The results in Fig. 7 are
with MIXED_HEIGHTS$=$’F’ . With MIXED_HEIGHTS$=$’T’ one obtains the results of
Fig. 24. The results without mixed heights and those from Basilisk are added
for comparison. There is a difference remaining with the Basilisk computations
than we have not yet been able to explain.

Figure 24: Maximum $L_{\infty}$ error norm in two dimensions for the curvature
estimated for a cylinder using the height function method in Paris and
Basilisk. Two Paris results are shown, one with the mixed curvature option and
one without. Averaging is used in both cases. Using the mixed-cell option
yields less accurate results than not using it. Figure 25: Maximum
$L_{\infty}$ error norm in three dimensions for the curvature estimated for a
sphere using the height function method in Paris and Basilisk. Two Paris
results are shown, one with the mixed curvature option and one without.
Averaging is used in both cases. Using the mixed-cell option yields less
accurate results than not using it.

###  A.4 Averaging scheme

A new loop over all cells cut by the interface is started. If both schemes
above have failed in the current cell, an average is performed over
neighboring cells that have been successful by either method in the first
pass. For each cell $\Omega_{i,j,k}$, the cubic set of neighbors

$B_{i,j,k}=\\{\Omega_{l,m,n}|i-1\leq l\leq i+1,j-1\leq m\leq j+1,k-1\leq n\leq
k+1\\}$

is defined. If at least one of the cells in $B_{i,j,k}$ has been successful,
the resulting curvature in $\Omega_{i,j,k}$ is the average curvature of these
successful cells.

###  A.5 Second pass: centroid fit

A final loop on cells $\Omega_{i,j,k}$ is performed. If all previous
approaches have failed or are not set to be used, then one falls back to a
fitting of centroids. In each cell of the cubic set $B_{i,j,k}$ containing an
interface, the cell centroid is computed using the vof_functions microlibrary
included in the code (implementing the method in [28]). Except for very small
fragments, the interface must find at least five neighboring cells in addition
to the current cell. This gives six points with which to fit the parameters in
expression (51).

In some cases, the fit fails because the least-square linear operator is not
invertible. The code then reports in a statistical manner these failures and
flags the cell as having an uncomputable curvature. A zero surface tension
force is then added to the momentum.

###  A.6 Comparison with other implementations of height-function curvature

The accuracy contrast when modifying the mixed-heights option is even more
dramatic in 3D, see Fig. 25. There is a striking drop in the $L_{1}$ error
near 13 grid points per diameter. This drop is due to the averaging step in
A.4. Suppressing the averaging reverts the results to accuracies comparable to
those of Basilisk or slightly better.