A phase inversion benchmark

for multiscale multiphase flows

S. Vincent§,†,a, L. Osmar†, J.-L. Estivalezes♠, S. Zaleski+,
F. Augusteo,∇, W. Aniszewski∗,+, Y. Ling+,b, T. Ménard∗,
A. Pedronoo, J. Magnaudeto, J.-P. Caltagirone†,
A. Berlemont∗, S. Popinet+
<EMAIL_ADDRESS>(corresponding author)
§ Université Paris-Est Marne-La-Vallée and CNRS, Laboratoire Modélisation et
Simulation Multi Echelle (MSME), UMR 8208, F-77454, Marne-La-Vallée, France
† Université de Bordeaux and CNRS, Institut de Mécanique et Ingénierie (I2M),
UMR 5295, F-33607 Pessac, France
o Université de Toulouse and CNRS, Institut de Mécanique des Fluides de
Toulouse (IMFT), UMR 5502, F-31400 Toulouse, France
+ Sorbonne Université and CNRS, Institut Jean Le Rond d’Alembert, UMR 7190, F-75005 Paris, France
∗ Université de Rouen and CNRS, Complexe de Recherche Interprofessionnel en
Aérothermochimie (CORIA), UMR 6614, F-76801 Saint-Etienne-du-Rouvray Cedex,
France
b Baylor University, Department of Mechanical Engineering, Waco, TX 76706, USA
♠ ONERA, The French Aerospace Lab, F-31055 Toulouse, France
∇ CERFACS, Centre Européen de Recherche et de Formation Avancée en Calcul
Scientifique, F-31057 Toulouse, France

###### Abstract

A series of benchmarks based on the physical situation of “phase inversion”
between two incompressible liquids is presented. These benchmarks aim at
progressing toward the direct numerical simulation of two-phase flows. Several
CFD codes developed in French laboratories and using either Volume of Fluid or
Level Set interface tracking methods are utilized to provide physical
solutions of the benchmarks, convergence studies and code comparisons. Two
typical configurations are retained, with integral scale Reynolds numbers of
$1.37\,10^{4}$and $4.33\,10^{5}$, respectively. The physics of the problem are
probed through macroscopic quantities such as potential and kinetic energies,
interfacial area, enstrophy or volume ratio of the light fluid in the top part
of the cavity. In addition, scaling laws for the temporal decay of the kinetic
energy are derived to check the physical relevance of the simulations. Finally
the droplet size distribution is probed. Additional test problems are also
reported to estimate the influence of viscous effects in the vicinity of the
interface.

keywords: phase inversion, multiphase flows, benchmark, Volume Of fluid, Level
Set, Ghost Fluid

##  1 Introduction

The physical validation or verification of CFD tools devoted to multiphase
flows and non miscible fluids, (by non-miscible we mean that interfaces remain
sharp), is a major concern. A lot of work was carried out in the 90s to
validate interface tracking algorithms in cases where the velocity field is
given analytically [1] [2]. Comparisons have also been achieved with available
experiments or predictions of linear stability theories [3]. However, in
almost all of these situations, the topology of the interface is rather simple
and the flow problem lies in a single multiphase characteristics such as the
continuity of stresses across the interface with the two-phase Poiseuille flow
[4], the Laplace law [5], the oscillation modes of a free surface [6] [5], the
head-on coalescence of two drops [7] [8], the dam break problem [9] or the
break-up of a low-velocity round jet due to the Rayleigh-Plateau instability
[10] to mention just a few.

None of these problems puts into play the interaction between interfaces and
turbulence or unsteady flows which is a key point in many realistic cases
involving a two-phase dynamics. Indeed, such multiphase interactions occur in
many environmental and industrial applications. Let us just mention the
atomization of a liquid jet in an engine[11, 12, 13], phase separators in
chemical engineering processes, the boiling crisis in nuclear reactors[14],
complex bubble flow including the flow past Taylor bubbles in pipes[15],
droplet or asteroid impact[16], wave impact on structures or the ubiquitous
phenomenon of wave breaking at the ocean surface[17]. In most multiphase flows
of environmental or engineering relevance, the multiscale character of the
flow is a key issue to be handled by computational approaches. The interfacial
length scales are associated on the one hand to large interface structures
such as jets, films or large drops. On the other hand, small interfacial
scales also exist, corresponding to a small-scale dispersed phase, i.e. small
bubbles or droplets. These two widely distinct families of interfacial scales
interact in a nonlinear way through ligament break-up or drop/bubble
coalescence and the unsteady or turbulent nature of the carrier fluid motion
plays a crucial role in these interactions. The definition of a reference
benchmark for these multiscale sharp interface problems appears to be an
important issue.

Figure 1: Example of simulated phase inversion in a closed box [18]
illustrating the multiscale character of the interfacial flow.

Homogeneous Isotropic Turbulence (HIT) often serves as a reference problem in
the verification of single-phase DNS codes. In analogy, one would like to
define a reference problem for the verification of CFD multiphase flow codes
devoted to non miscible fluids. Our first idea was to try to define the
reference problem on the basis of a well documented experiment, such as
atomization of liquid sheets or jets. However, the quality of modeling and
simulation of atomization is very sensitive to inlet forcing conditions, in
particular due to upstream turbulence, boundary layers or inlet pump
oscillations [19, 20] . Although numerous experimental studies can be found in
the literature [9, 21, 22] they are plagued by two difficulties that make them
unsuitable as references cases: either they have too simple flow structures
(as for dam break on a dry bottom) or they are not sufficiently well
documented or instrumented to provide all the information required for
defining boundary and initial conditions for simulations. Moreover many
experiments lack instantaneous pressure, velocity and interface measurements
(dam break on a wet bottom [23]). As a consequence, we have decided to define
a synthetic reference case unrelated to a specific experiment, although an
experiment in this configuration could in principle be considered (the setup
resembles that of an oil and water phase separator or settling pond). In this
reference case initial and boundary conditions are simple and well known. At
the same time, we wish to have complex flow features (turbulence, interface
rupture, coalescence, multi-scale flow characteristics). A setup we believe
satisfies all these requirements is a cubic outer geometry and a simple
arrangement of the phases both as initial condition and final equilibrium
state of the dynamics, so that grid quality and its influence on the Navier-
Stokes solver performance will not be a problem. At the same time, there are
advantages to consider sufficiently large Reynolds number ${\rm Re}$ and Weber
number ${\rm We}$ to have the multiphase equivalent of “fully developped
turbulence”, a state in which the statistical features of the flow on large
and intermediate scales become relatively independent of ${\rm Re}$ and ${\rm
We}$ as those numbers further increase. Such fully developped multiphase
turbulence is the relevant flow state with respect to common natural sciences
or engineering situations.

The simple cubic geometry idea leads naturally to the phase inversion of a
cube of light fluid inside a larger cube of a heavier fluid contained in a
closed box. As desired, with high enough ${\rm Re}$ and ${\rm We}$ this flow
configuration is typical of fully developped turbulent multiscale flows in
which a broad range of interfacial scales (i.e. those characterizing the
typical size of the interfaces separating the two fluids) coexist: as a result
of multiple break-up and coalescence events, large-scale sloshing motions are
observed together with small dispersed droplets interacting with larger drops.
An illustration is given in Figure 1 taken from [18], [24]. As described in
Figure 2, the initial condition is geometrically straightforward and the
solution at large times is known, i.e. in the final stage the two fluids are
superimposed in a stable manner (with all the light fluid in the top part of
the box) and are separated by a single horizontal interface. A parametric
study in which various physical parameters are varied can easily be performed
by changing the density or viscosity ratios as well as the surface tension.
From a numerical point of view, the physical parameters can be defined in such
a way that a numerical convergence study may be achieved until a true direct
numerical simulation (DNS) is obtained, i.e. all the flow and interfacial time
and space scales are resolved. The sensitivity of the results with respect to
the grid density or to the numerical methods can also be examined. Laminar as
well as turbulent conditions can be generated. Regimes corresponding to either
thinly dispersed flows or to large-drop topologies can be selected by properly
choosing the physical characteristics of the two fluids.

Among the wide number of references existing on the subject of multiphase flow
modeling and simulation, it appears to us that there are two categories of
work, those concerning mixing and diffuse evolution of initially discontinuous
distribution of fluids and those associated to interfaces remaining sharp
during time, i.e. immiscible fluids. The first category has investigated for
example the interaction between Rayleigh-Taylor (RT) instabilities and
turbulence in numerous works [25, 26, 27, 28] and a benchmark has been
published in 2004 by Dimonte and co-workers [29]. This reference work involves
physical phenomena similar to those occurring in the phase inversion
configuration proposed here. Much attention has been paid to characterizing,
with various models and numerical methods, the effects of initial
perturbations of the interface, volume fraction profiles, self-similar bubble
dynamics or mixing behavior. However, the major difference is related to the
immiscible nature of fluids and discontinuous character of interface (on a
macroscopic point of view). It has been clearly demonstrated [28] that the
diffuse character of the fluids has a major effect on the global dynamics of
the turbulent two-phase flow in terms of spike position for example and that
this discrepancy was increasing during time. In diffuse modeling of RT flows,
simulations mainly involve gas stratifications whereas our focus here is on
immiscible fluids (with at at least to one liquid phase in the considered
problem). Our benchmark thus belongs to the second category of problems for
which RT simulations also exist [30, 31]. Intermediate modeling between
diffuse and sharp interfaces have also been proposed for simulating Rayleigh-
Taylor instabilities [32]. The major concern was the model and numerical
methods and there was no consideration of the turbulence characteristics of
the flows. To our knowledge, neither benchmark nor detailed simulation works
have been published about RT interaction with turbulence in the framework of
sharp liquid interfaces and immiscible fluids. In that sense, the present work
is a relevant and new contribution for the CFD community of multiphase flows
with sharp interfaces.

The convergence study is the first original outcome of this work. Indeed,
neither detailed nor accurate experimental measurements seem to exist for the
phase inversion problem considered here. Therefore it is of fundamental
interest to validate the quality of the various technical s.png of the
computations, namely the numerical schemes, the Navier-Stokes solvers, etc.
Our second objective is to determine whether or not the interface tracking
techniques and the numerical methods required to discretize the two-phase
Navier-Stokes equations (e.g. the definition of the two-phase viscosity as a
function of those of the two fluids) can provide consistent solutions for the
multiphase and multiscale flow problem at hand. To meet both objectives, a
comparison of the performance of five distinct in-house DNS multiphase flow
codes on the reference or benchmark case is provided in what follows. It
highlights the capabilities but also the limitations of up-to-date numerical
methods and interface models to tackle multiscale, multiphase, flow problems.

##  2 Definition of the phase inversion benchmark

The phase inversion configuration has previously been proposed or used [33]
without the objective of analyzing the physical aspects of the multiscale
character of the flow, nor comparing various interface tracking methods. A
benchmark based on this flow situation was also recently performed using VOF
methods [31]. The simulation of this configuration was used to obtain a priori
estimates of subgrid terms that appear and must be modeled when LES is used in
the framework of two-phase flows of immiscible fluids [34]. A new selection of
the physical parameters is achieved here (see Tables 1,2 and 3). With the very
high Re and La numbers achieved (Table 3) viscosity is considerably smaller
than what would be needed to obtain fully-resolved simulations in multiphase
as well as in single phase flow. (Note that other computations were performed
in [35] [33] [18] [24].) The two cases that are presented in this work are a
good illustration of an atomization process at large Reynolds number. We
cannot hope to perform a Direct Numerical Simulation in the sense of a
simulation with all length scales resolved. Clearly, our simulations without
explicit turbulent or subgrid droplet modeling, are finely-resolved implicit
Large Eddy Simulations.

###  2.1 Geometry, initial and boundary conditions

Figure 2: Sketch of the phase inversion problem in a closed box.

As described in Figure 2, an initial cubic blob of light liquid, referred to
as fluid 1, is placed in the bottom part of a cubic box filled with a heavier
liquid, referred to as fluid 2. The size of the box is $(H,H,H)$, while the
size of the blob of light fluid is $(H/2,H/2,H/2)$. All outer walls are
considered as free-slip impermeable walls, so that the normal velocity is zero
and the tangential components obey a symmetry condition. A $\pi/2$ static
contact angle is implicitly assumed on the walls. Gravitational acceleration
is chosen as ${\bf g}=(0,-9.81,0)$. Besides the simplicity of the initial and
final configurations of the phase inversion problem, another of its advantages
is that it provides the possibility to observe multiple coalescence and break-
up events, although that introduces an additional complexity since in a
simulation coalescence is conditioned by numerical rather than physical
properties. In this setup we may conveniently  observe fragmentation of
ligaments in a kind of atomization process, while not depending on complex
initial or forcing conditions such as those encountered in real atomization
configurations [36], [37]. The main disadvantage of this benchmark case is
that no experimental data are available.

###  2.2 Physical parameters of the test cases

Characteristics of the fluids are varied in order to investigate various
multiscale interface topologies. The fluid properties control these
conditions, which may be expressed in terms of Re, and We. In order to define
these numbers we need to define a velocity scale. We obtain it from energy
considerations. In all of what follows, we use the characteristic or color
function $C$ that is defined with respect to fluid 1, so that $C=1$ in fluid 1
and $C=0$ in fluid 2. The potential energy of phase $n$ is

$E_{p,n}=\int_{\Omega}C_{n}\rho_{n}gydV$ (1)

where $C_{k}$ is the characteristic function of phase $n$: $C_{1}=C$,
$C_{2}=1-C$ and with $g=\|{\bf g}\|$). The kinetic energy of phase $n$ is

$E_{k,n}=\frac{1}{2}\int_{\Omega}C_{n}\rho_{n}{\bf u}^{2}dV$ (2)

Due to the simple topology of the interface in the initial and final s.png the
potential energies are trivially found. At $t=0$,
$E^{i}_{p,1}=\rho_{1}g{H^{4}}/{32}$ and $E^{i}_{p,2}=15\rho_{2}g{H^{4}}/{32}$
while for $t\rightarrow\infty$, $E^{f}_{p,1}=15\rho_{1}g{H^{4}}/{128}$ and
$E^{f}_{p,2}=49\rho_{2}g{H^{4}}/{128}$. An upper bound on the kinetic energy
is

$E_{k,1}+E_{k,2}\leq
E^{i}_{p,1}+E^{i}_{p,2}-E^{f}_{p,1}-E^{f}_{p,2}=\frac{11}{128}(\rho_{2}-\rho_{1}){gH^{4}}$
(3)

and an upper bound on the velocity is $||{\bf
u}||_{2}<(2E_{k}/\rho_{1})^{1/2}$. Identifying the velocity with the
upperbound we obtain the velocity scale

$U_{K}=\frac{\sqrt{11}}{8}\sqrt{\frac{\rho_{2}-\rho_{1}}{\rho_{1}}gH}$

The dimensionless numbers are defined based on $U_{K}$ and properties of fluid
2,

$\displaystyle{\rm Re}=\frac{\rho_{2}HU_{K}}{\mu_{2}}\,,$ (4)
$\displaystyle{\rm We}=\frac{\rho_{2}HU_{K}^{2}}{\sigma}\,,$ (5)

As is customary in atomization problems, the Reynolds and Weber numbers may be
combined in a way that eliminates the velocity scale, yielding the Laplace
number La=$\sigma\rho_{2}H/\mu_{2}^{2}={\rm Re}^{2}/{\rm We}$.

To allow comparison with previous work on this setup we also define the same
velocity scale as in previous studies by some of us [18, 24, 34]
$\displaystyle{U_{g}}={\frac{\rho_{2}-\rho_{1}}{\rho_{1}}\sqrt{\frac{Hg}{2}}}$
and a characteristic time scale $t_{c}=H/(2U_{g})$.

The fluid properties correspond approximately to oil (fluid 1) and water
(fluid 2). These configurations are interesting in so far as they provide
interface and turbulent characteristics similar to those of jet atomization
with simpler and more easily controlled initial conditions, and are as a
result more convenient for numerical simulations. These liquid-liquid two-
phase flows make it possible to select easily a range of Reynolds numbers
smaller, or even equivalent, to those encountered in real atomization
processes while exhibiting equivalent dynamics. The characteristics of fluid 1
and related dimensionless numbers are given in Tables 1,2 and 3. Case 1
corresponds to a configuration with moderate fragmentation. Case 2 has
increased size and thus larger Reynolds number (by a factor of 30 compared to
case 1). The Weber number is also higher in case 2, (by a factor of 10
compared to case 1) but does not grow as rapidly as the Reynolds number since
surface tension was also artificially increased. This was done specifically to
avoid having an exceedingly large Weber number which would lead to excessive
fragmentation of the interface. Nevertheless both cases involve large inertial
effects and low capillary forces resulting in the generation of a large number
of droplets.

$\rho_{1}$ | $\rho_{2}$ | $\mu_{1}$ | $\mu_{2}$ | $g$
---|---|---|---|---
($kg.m^{-3}$) | ($kg.m^{-3}$) | ($Pa.s$) | ($Pa.s$) | $(m.s^{-2})$
900 | $1000$ | 0.1 | $0.001$ | 9.81
Table 1: Common fluid and physical properties Case | $\sigma$ | $H$
---|---|---
| ($kg.s^{-2}$) | $(m)$
1 | 0.045 | 0.1
2 | 0.45 | 1
Table 2: Fluid and physical properties that depend on the case considered. Case | ${\rm Re}$ | ${\rm We}$ | La=${\rm Re}^{2}/{\rm We}$
---|---|---|---
1 | $1.37\,10^{4}$ | $41.6$ | $4.5\,10^{6}$
2 | $4.33\,10^{5}$ | 416 | $4.5\,10^{8}$
Table 3: Dimensionless numbers for the cases considered.

###  2.3 Macroscopic quantities of interest

Several macroscopic quantities characterizing the evolution of the flow field
are of primary interest. Among these physically relevant quantities we select
the following.

  * •

A measure of the amount of fluid 1 that has reached the upper part of the box.
The upper part is the box $\Omega_{up}$ of dimension $(H,H,H_{2})$ at the top
of domain. As $t\rightarrow\infty$ it is the exact location of fluid 1. Given
the initial volume of fluid 1, we have $H_{2}={H}/{8}$ and the ratio

$\displaystyle{R_{2}=\frac{1}{V_{2}}\int_{\Omega_{up}}C_{2}dV}$ (6)

is a relevant quantitative representation of the phase separation process
which has the asymptotic value $\displaystyle{R_{2}=1}$. In the benchmark,
values of (6) will be saved at each time step of the simulations. Their
evolution will help us estimate the typical time of phase separation time and
characterize the macroscopic dynamics of the problem.

  * •

The potential and kinetic energies will be stored at each time step. They will
help monitor the conversion of potential energy into kinetic energy and
provides a characterization of the density stratification.

  * •

The time evolution of the volume integral of the enstrophy in both fluids will
also be recorded. This quantity is defined as $E_{r,n}=\int_{\Omega}C_{n}{\bf
w}^{2}dV/2$, with ${\bf w}=\nabla\times{\bf u}$ denoting the vorticity.

  * •

The time evolution of the interfacial area will also be recorded, being a
sensitive marker of grid convergence. This parameter is also relevant in
problems involving interfacial energy or mass transfer [38], [39].To estimate
it, we introduce a binary presence function $C^{b}$ such that $C^{b}=0$ if
$C<0.5$ and $C^{b}=1$ elsewhere. This function can also be easily built from
the level set function. The interfacial area $\Sigma_{int}$ is finally
approximated as

$\Sigma_{int}\approx\sum_{\\{i,j,k\colon\|\nabla C^{b}_{i,j,k}\|\neq
0\\}}\Delta x\Delta y\,.$ (7)

##  3 Model and numerical methods

As is now well established, incompressible two-phase flows involving fluid-
fluid interfaces and Newtonian fluids can be modeled by a single set of
incompressible Navier-Stokes equations with phase-specific density and
viscosity and possibly extra interfacial forces (e.g. the capillary force),
together with the transport equation of the phase function $C$. The resulting
model takes implicitly into account the mass and momentum jump relations at
the interface [40] [41], whereas the continuity of the fluid-fluid and fluid-
solid interfaces are taken into account by the $C$ equation. The entire set of
equations reads:

$\displaystyle\nabla\cdot{\bf u}$ $\displaystyle=$ $\displaystyle 0\,,$ (8)
$\displaystyle\rho\left[\frac{\partial{\bf u}}{\partial t}+({\bf
u}\cdot{\nabla}){\bf u}\right]$ $\displaystyle=$ $\displaystyle-\nabla
p+\rho{\bf g}+\nabla\cdot\left[\mu\left(\nabla{\bf u}+\nabla{\bf
u}^{T}\right)\right]+{\bf F}_{st}\,,$ (9) $\displaystyle\frac{\partial
C}{\partial t}+{\bf u}\cdot\nabla C$ $\displaystyle=$ $\displaystyle 0\,,$
(10)

where $p$ is the pressure, ${\bf F}_{st}$ is the interfacial force per unit
volume and $\rho$ and $\mu$ are the local density and viscosity of the two-
phase medium, respectively.

Capillary effects are inserted in the source term ${\bf F}_{st}$ in the form
${\bf F}_{st}=\sigma\kappa{\bf n}_{i}\delta_{i}$, where $\sigma$ denotes the
surface tension, $\kappa$ is the local mean curvature of the interface, ${\bf
n}_{i}$ is the unit vector normal to the interface and $\delta_{i}$ is the
interface Dirac function [42].

The above one-fluid model is the classical model for multiphase incompressible
flow with sharp interfaces and surface tension. Localizing the interface
requires solving the additional transport equation for the characteristic
function $C$. This can be performed using the volume fraction (also noted $C$
in Volume-of-Fluid methods) or using a continuous level-set function $\phi$
such that $\phi=0$ at the interface and $\phi>0$ (resp. $<0$) in fluid 2
(resp. 1). In this case, the characteristic function is obtained as
$C=H(\phi)$ where $H$ is the Heaviside function[43].

Note that the above equations are clearly distinct from those appearing in
various models of single or multiphase turbulence, that may involve Reynolds
stresses, terms for unresolved eddies or turbulent viscosities and
diffusivities, and non-sharp interfaces.

###  3.1 The Thétis code

Thétis is a CFD code developed in the TREFLE Department of the I2M laboratory.
It solves the one-fluid Navier-Stokes equations discretized with implicit
finite-volumes on an irregular staggered Cartesian grid. A second-order
centered scheme is used to approximate the spatial derivatives while a second-
order Euler or Gear scheme is used for the time integration [44]. All terms
are written at time $(n+1)\Delta t$, except the inertial term which is
expressed in the following semi-implicit manner:

${\bf u}^{n+1}\cdot\nabla{\bf u}^{n+1}\approx\left(2{\bf u}^{n}-{\bf
u}^{n-1}\right)\cdot\nabla{\bf u}^{n+1}.$ (11)

It has been shown that this approximation allows to reach second-order
convergence in time [45]. The coupling between velocity and pressure is
ensured by using an implicit algebraic adaptive augmented Lagrangian method
[46]. The augmented Lagrangian methods used in this work are independent on
the chosen discretization and could for instance be implemented in a finite-
element framework [47]. In two dimensions, the standard augmented Lagrangian
approach [48] can be used to deal with two-phase flows as direct solvers [49]
are efficient in this case. However, as soon as three-dimensional problems are
considered, the linear system resulting from the discretization of the
augmented Lagrangian terms has to be treated with a BiCG-Stab II solver,
preconditioned by a Modified and Incomplete LU method [50]. As for the
interface tracking and advection equation of $C$, two different volume of
fluid (VOF) methods have been implemented in Thétis [51] [52]. They are
evaluated here. The above numerical methods and the one-fluid model have been
validated in previous works, e.g. [53] and [54].

###  3.2 The Gerris Flow Solver code

The Gerris Flow Solver (GFS) is a free code implementing finite volume solvers
on an octree adaptive grid together with a piecewise linear VOF interface-
tracking method. In this work, the Navier-Stokes equations for an
incompressible two-phase flow with constant surface tension are solved. The
simulations reported in this paper use a VOF method of the piecewise linear
type, in which the interface segments are reconstructed using the mixed-
Youngs-centered (MYC) approximation [55]. While not the most accurate for very
fine grids, the MYC approximation is easily implemented when using only
information from the nearest-neighbour cells, an important advantage when
domain-decomposition on octrees is used. Advection of $C$ is performed using
the Lagrangian-Explicit method first published by Li [56] and discussed in
[57, 58, 59]. While not volume-conserving to machine accuracy, it has very
good volume and mass conservation properties. Surface tension is a vexing
question in multiphase flow. As density and viscosity ratios become very large
or very small, the simulations become increasingly difficult with standard
methods [59]. The problem has been considerably improved in Gerris as a result
of the use of Height-Function methods [5] and a so-called balanced-force
algorithm [60, 61]. For a review of surface-tension methods including a
discussion of the differences between those used in Gerris and in other codes,
see [62]. Gerris has been used with success in two-dimensional [63] and three-
dimensional [64, 65] atomization and droplet impact studies.

###  3.3 The Jadim code

JADIM is a versatile code developed for a number of years at IMFT (Institut de
Mécanique des Fluides de Toulouse). In JADIM, the time evolution of two-phase
or three-phase flow is obtained using the one-fluid formulation of the Navier-
Stokes equations. The momentum equations are discretized on a staggered
orthogonal grid using a finite-volume approach. Spatial discretization is
performed using second-order centered differences. Time-advancement is
achieved through a third-order Runge-Kutta algorithm for advection/source
terms and a Crank-Nicolson algorithm for viscous stresses. Incompressibility
is satisfied at the end of each time step through a projection method. More
details may be found in [38]. The resulting code is second-order accurate in
both time and space for single-phase flows.

In two- and three-phase configurations, a VOF method with no interface
reconstruction is used. The advection equation for $C$ is solved using a flux-
corrected scheme split into successive one-dimensional s.png [66]. Owing to
the splitting procedure, the overall transport scheme is not rigorously
conservative. Therefore a local mass error control which improves upon the
global control technique described in [56] is employed. The corresponding
strategy is based on the detection of drop or bubbles (through an index
function and a topological monitoring of each bubble or drop) and on an
iterative solution which modifies the volume fraction within the transition
regions (i.e. those in which the volume fraction is neither zero nor one) in
order to keep the volume of each bubble or drop constant. Since no interface
reconstruction step is involved, smearing of interfaces frequently occurs in
high-shear regions. A specific strategy is employed to keep this smearing
within reasonable bounds. For this purpose, the velocity at nodes crossed by
the interface is modified to keep the thickness of the interfacial region
constant. Details may be found in [4, 66]. In the most difficult cases (e.g.
break-up and coalescence), this approach is supplemented by an antidiffusion
technique in which the local volume is redistributed in the direction normal
to the interface so as to eliminate spurious values of the volume fraction
outside interfacial regions.

The capillary force is represented using a modified version of the Continuum
Surface Force (CSF) model [42]. The modification used in JADIM consists in an
original manner of evaluating the interfacial area. It allows us to define the
interfacial region as that in which the capillary force takes non-zero values
(not as that in which the volume fraction takes intermediate values) and
improves the control of the thickness of the transition region.

###  3.4 The Archer code

Archer is a CFD code developed in the CORIA Laboratory mainly devoted to
multiphase flows. It has been previously applied to compute atomization [36],
vaporization and mixing [67] and other complex interfacial flows [31]. The
Level Set (LS) method [43] is used for tracking interfaces and shapes. It is
based on a continuous distance function $\phi$ defined as the signed distance
between any point of the domain and the interface. Similar to the volume
fraction $C$, $\phi$ obeys a pure advection equation.

To avoid singularities in the $\phi$ field, the fifth order conservative WENO
[68] scheme is applied to discretize convective terms. When the LS advection
is carried out, high velocity gradients can cause wide spreading or stretching
of $\phi$ which then no longer remains a distance function. A redistancing
algorithm [69] is thus applied at every time step to restore the distance
property of $\phi$, i.e. $|\nabla\phi|=1$. Advancement of the $\phi$-equation
and the redistancing algorithm can induce mass loss in under$-$resolved
regions [31]. This is the main drawback of Level Set methods. To improve mass
conservation, the Coupled Level-Set Volume of Fluid (CLSVOF) method [70] is
used. The main idea of this method is to benefit from the advantages of each
tracking strategy: minimize the mass loss using VOF and keep a fine
description of interface properties with the smooth LS function. (A study of
the performance of various VOF variants was recently published [31] using
Archer in the CLSVOF context; the phase inversion problem is also discussed
therein.) The coupling between LS and VOF is maintained using a correction
scheme based on geometric reconstruction of both VOF and LS interfaces [34].
The LS method itself is coupled with a projection method for the
incompressible Navier$-$Stokes equations, where the density and the viscosity
depend on the sign of the LS function (with appropriate interpolations used in
interfacial cells). To finalize the description of the two-phase flow, jump
conditions across the interface are taken into account using the Ghost Fluid
(GF) method. In the GF approach, ghost cells are defined on each side of the
interface [71] [72] and appropriate schemes are applied to the jump of each
quantity. As defined above, the interface is characterized through the
distance function, and jump conditions are extrapolated on some nodes on each
side of the interface. Following the jump conditions, the discontinued
functions are extended continuously and then derivatives are estimated.
Heaviside functions are designed according to the distance function in order
to provide a characteristic function for the two-phase medium. They follow the
functions proposed by [73] according to formula (6) in its publication.
Concerning discrete Dirac functions, they are built in a similar way and
regularized by a properly set cosine function. These techniques have been
presented and validated in previous works, e.g. [36].

###  3.5 The DyJeAT code

DyJeAT (Dynamics of Jet ATomization) is an in-house computational fluid
dynamics library developed at ONERA. Velocity/pressure coupling is ensured
with classical projection methods [74] [75]. In the same way as the ARCHER
code, the Ghost Fluid method [76] is used to deal with surface tension forces
as well as density and viscosity jumps. The Level-Set (LS) approach is also
used for tracking interfaces and LS-VOF coupling is also used to improve mass
conservation. An adaptive mesh refinement (AMR) capability has been added to
DyJeAT in order to concentrate computational effort where it is most needed.
It is based on block-structured AMR. Details of implementation and
verification can be found in [77]. As the physics of atomization processes
involve many spatial scales, generating a wide range of droplet sizes with
large density and viscosity ratios, each scale has to be resolved with an
appropriate method to ensure the conservation of mass and momentum and satisfy
the jump conditions across the interface. To address these problems, the LS-
VOF approach is coupled to a particle tracking technique ($via$ an improved
Eulerian-Lagrangian coupling) to capture the droplet dynamics [78].

###  3.6 Pros and cons of different codes

Code | Interface | Navier-Stokes | Capillary | Grids | Solvers
---|---|---|---|---|---
| tracking | velocity-pressure | forces |  |
|  | coupling |  |  |
Archer | CLSVOF | Projection | Ghost | Cartesian | BiCGStab
|  | method | fluid | staggered | Multigrid
DyJeat | Level set | Projection | Ghost | Cartesian | BiCGStab
|  | method | fluid | staggered | Multigrid
Gerris | Geometrical | Projection | CSF with | Collocated | BiCGStab
| VOF | method | height | and AMR | Multigrid
|  |  | functions |  |
Jadim | Implicit | Projection | CSF | Cartesian | BiCGStab
| VOF | method |  | staggered | Multigrid
Thetis | Geometrical | Augmented | CSF with | Cartesian | BiCGStab
| VOF | Lagrangian | smooth VOF | Staggered | ILU
Table 4: Summary of the numerical methods used in the various codes.

All the codes used in the present work are based on finite volumes but with
different variants of the methods. Obviously, each code can lead to a
different result for a given grid, due to variations in the methods used.
Table 4 is given to clarify the comparisons of the model and numerical method
used. Schematically, the codes cover the broad range of numerical methods that
can be found in literature for simulating multiphase flows using the sharp
interface approximation described above. This approach is sometimes called
direct or detailed to contrast it with other approaches that depart from the
sharp interface approximation. Four of the codes use staggered grids while one
(Gerris) is using a collocated grid. Staggered grids have several advantages:
a more accurate pressure solution, and the avoidance of spurious “red-black”
velocity oscillations. However Gerris has demonstrated high accuracy on test
cases [5, 65] and ten years of experience with Gerris show almost no
occurrence of spurious oscillations. Four codes are using projection velocity-
pressure coupling that contains a time splitting error while one code (Thetis)
uses an exact augmented Lagrangian approach. This technique is exact if the
residual of the iterative solver is zero at machine accuracy. In terms of
interface tracking, Gerris and Thetis are using a geometrical VOF-PLIC
technique while Jadim utilizes a FCT scheme to directly approximate the
hyperbolic advection equation on the VOF. DyJeat uses the Level Set method
whereas Archer couples this method with VOF-PLIC in order to improve mass
conservation of Level Set and also better inertia treatment in the momentum
conservation equation. Concerning surface tension, the most accurate
approaches are the Ghost fluid that uses jump relations at the interface and
the CSF with height functions [5]. Thetis and Jadim uses less sophisticated
capillary force approximations that are known to generate more spurious
currents than other methods. To finish with numerical methods, all projection
techniques are solved with iterative BiCGStab II solver preconditionned with a
multigrid algorithm while Thetis and the coupled augmented Lagrangian method
uses BiCGStab II and an incomplete LU preconditionner, because the linear
system is not symmetric and more difficult to solve than with the projection.
The major drawback of the Thetis solver is that its parallelization is keeping
a good speed-up until $1000$ processors but it can hardly handle in the
present form more processors due to the ILU preconditionner. On the contrary,
the multigrid preconditionner used in projection is nicely extendable in
parallel until $100000$ processors.

##  4 Benchmark of phase inversion

The teams that have participated to the present benchmark have each run their
own code. All the codes were run on $256^{3}$ grids. Moreover one of the codes
(DyJeAT) was run on a sequence of much finer grids in order to provide a test
of convergence. Cross-comparisons lead us to conclusions on the physical
behavior of the numerical methods and also on the physical meaning of the
simulations.

###  4.1 Case 1: a phase inversion problem with few fragmentation events

Parameter | Value | Units
---|---|---
$t^{*}={t}/{t_{c}}$ | $\displaystyle{\frac{t}{0.643}}$ | -
$E_{p,1}^{f}$ the potential energy in fluid 1 for $t\rightarrow\infty$ | $0.1035$ | J
$E_{p,2}^{f}$ the potential energy in fluid 2 for $t\rightarrow\infty$ | $0.3755$ | J
$\displaystyle{E_{k,1}^{*}=\frac{E_{k,1}}{(1/16)\rho_{1}U_{g}^{2}H^{3}}}$ | $\displaystyle{\frac{E_{k,1}}{0.000341}}$ | -
$\displaystyle{E_{k,2}^{*}=\frac{E_{k,2}}{(1/16)\rho_{2}U_{g}^{2}H^{3}}}$ | $\displaystyle{\frac{E_{k,2}}{0.000378}}$ | -
${{H^{3}}/{8}}$ the final volume of fluid 1 in the top part of the box | $0.000125$ | $m^{3}$
Maximum of enstrophy in fluid 1 (DyJeAT on a $512^{3}$ grid) | $0.0733$ | $m^{3}.s^{-2}$
Maximum of enstrophy in fluid 2 (DyJeAT on a $512^{3}$ grid) | $1.3759$ | $m^{3}.s^{-2}$
${{3H^{2}}/{4}}$ the initial interfacial area of the drop of fluid 1 | $0.0075$ | $m^{2}$
Table 5: Parameters used to define the dimensionless variables in case 1.

Figure 3: Grid convergence of kinetic energies for case 1 with DyJeAT.

Figure 4: Potential energies for case 1.

Figure 5: Kinetic energies for case 1 in log-linear (top) and log-log (bottom)
coordinates.

In this section we present the numerical results pertaining to case 1. The
characteristics of this case are summarized in Tables 1-3. The macroscopic
quantities defined in section 2.3 were computed by the five codes described in
the previous section. Dimensionless data are considered according to the
definitions indicated in Table 5.

The potential and kinetic energies are mostly dependent on the large scales of
the flow and can thus be dubbed “macroscopic” quantities. A convergence study
for these quantities was carried out on $128^{3}$, $256^{3}$ and $512^{3}$
grids with the DyJeAT code. The corresponding results are reported in Figure
3. There is still some difference between the $256^{3}$ grid and the $512^{3}$
grid for fluid 2 around time 2.5, although not in fluid 1. For the sake of
simplicity in what follows we will consider it to be the reference solution,
although this may not be a fully converged solution after time 2.5 . We note
that since as shown in Table 5 the kinetic energy is scaled by
${(1/16)\rho_{n}U_{g}^{2}H^{3}}$ for fluid $n$, the maximum kinetic of fluid 2
in the scaled variables is $8(u_{K}/u_{G})^{2}(\rho_{1}/\rho_{2})=11.1375$.
This is close to the maximum of the fluid 2 curve on Figure 3. The fact that
the kinetic energy is close to its upper bound also means that $U_{K}$ is
close to the $L_{2}$ norm of the velocity, and that our Reynolds and Weber
number estimates are accurate with the $L_{2}$ norm velocity.

The evolution of the potential and kinetic energies of the other codes is
plotted in Figure 4 and Figure 5, respectively, together with the reference
solution. Computations were run on a $256^{3}$ grid. The results provided by
all five codes on various grids are in approximate agreement with each other
and with the reference solution for these quantities. Around time $t^{*}=3$ we
see differences between all the codes during the first oscillation arc of the
potential energy. At later times, the oscillations of all the codes are
superposed, except for the two Thetis codes, and to a lesser degree the Jadim
code. For the kinetic energy the behavior is more complex as we see
differences between the behavior in fluid 1 and 2. However for both fluids, at
large times the kinetic energy predicted by the Archer code is the closest to
the reference solution. This may be due to the non-converged character of the
DyJeAT reference solution for $t^{*}>2.5$, which would make agreement with
Archer easier since both codes are at least partially Level-Set based. In
other words, Archer being the code most similar to DyJeAT will be the most
likely to agree with the reference. Figure 4 indicates that phase separation
is achieved after a dimensionless time equal to $15$. The dimensionless
frequency $f$ of the waves observed at the surface of the light fluid is about
$0.6$.

According to Figure 5, the kinetic energy decreases as $e^{-\alpha t^{*}}$ in
fluid 1, which is referred to as Stokes’ decay law (18), whereas it decays as
$t^{*-2}$ in fluid 2, which we refer to as the turbulent decay law (23). The
origin of these decay laws and the physical assumptions on which they are
based are discussed in Appendix A.

Figure 6: Evolution of the volume ratio of fluid 1 in the top part of the box
in case 1.

To characterize the phase separation on a macroscopic point of view, the
volume ratio of fluid 1 occupying the top $H/8$ part of the box is plotted in
Figure 6. A convergence study was carried out on $128^{3}$, $256^{3}$ and
$512^{3}$ grids with the DyJeAT code. This time the difference between the
$256^{3}$ grid and the $512^{3}$ grid is much smaller than the difference
between the $256^{3}$ grid and the $128^{3}$ grid. Again, it is observed here
that all codes predict evolutions close to each other and to the reference
solution. The code-to-code comparison indicates some phase shifts and various
damping rates at large times. The largest difference is found for the Thetis
code, followed by the Jadim code. The Archer code is even closer to the
reference and the Gerris code is closest. One possible explanation could be
the quality of mass conservation in the codes, since it would directly
influence the amount of fluid in the top part. Another explanation would be
the accuracy of energy conservation since it drives the oscillations seen. A
third explanation for the different behavior of the codes is their more or
less diffusive character. Level sets have the smoother indicator function
while VOF have the sharper, and Jadim has an intermediate behavior. Clearly
the ranking of the codes cannot be exclusively alloted to the conservation of
mass since all the VOF codes are supposed to do well in this respect but do
not perform equally well in the comparison. In addition, the effective initial
volume of light fluid is not initially strictly equal to the theoretical
volume of $H^{3}/8$. Also in some cases one observes liquid drops trapped on
the domain walls. Both of these effects explain the fact that numerically
obtained volume ratios tend to a value slightly less than 1.

These plots suggest that phase separation is almost completed after more than
$20$ time units, which is slightly longer than what was observed by examining
the potential energy.

Figure 7: Enstrophy for case 1. Figure 8: Enstrophy for case 1 in log-log
coordinates.

Figure 9: Convergence of enstrophy for case 1. Figure 10: Vorticity structures
and interface shape for case 1 - top: iso-60 vorticity magnitude (blue) and
iso $C=0.5$ interface (orange surface) - bottom: isoline of vorticity
magnitude and interface shape (black line) obtained with Thétis on two
different grids.

Figures 7 and 8 display the instantaneous values of the enstrophy in both
fluids. Although all codes were shown to agree well on other macroscopic
quantities, the magnitude of the enstrophy peak differs dramatically for each
of them. The location of this peak is found to occur at $t^{*}\approx 2$
(resp. $3$) for all codes in fluid 1 (resp. fluid 2). However, differences in
magnitude up to $80\%$ are observed, depending on the grid and the code. The
spatial convergence study of enstrophy carried out with DyJeAT is illustrated
in Figure 9. Although the enstrophy seems to converge in fluid 1, its peak
increases with the grid resolution in fluid 2 and no convergence is reached
even on the $512^{3}$ grid. Therefore, although primary moments such as
potential and kinetic energies or volume ratio of the light fluid suggest that
a real DNS was achieved, the analysis of higher-order moments such as
enstrophy invalidates this hope. Despite the large discrepancies noticed
between all codes, it is interesting to observe that they all provide a
$t^{*-3}$ temporal decay law which may be shown to correspond to a turbulent
scaling (see Appendix A).

The interface shape and isovalues of vorticity magnitude at $t^{*}=3$ (which
corresponds to the peak within fluid 2 in Figures 7 and 8) are reported in
Fig. 10. The maximum magnitude of the vorticity is found to be concentrated in
the near-wall regions and in regions where the large blob of light fluid that
has risen during the phase separation process undergoes high deformations.
Finer vortical structures are clearly captured by the finer grid inside fluid
2 (near the bottom corners for instance). When the large blob of light fluid
moves from one side of the box to the other, heavy fluid is trapped between
some of them and the wall, generating high shear (hence high vorticity)
regions. Such boundary layers are clearly not fully resolved with the grids
used in the present work, which explains the lack of convergence observed for
the volume-averaged enstrophy.
In single-phase turbulent flows, Direct Numerical Simulation is known to
require a number of grid points of the order of $Re^{9/4}$ [79] [80]. Provided
no scale smaller than the Kolmogorov microscale exists in the phase inversion
problem investigated here, and considering the integral scale Reynolds number
based on the parameters introduced in section 2.3, at least
$13700^{9/4}=2.03\cdot 10^{9}$ grid points would be necessary to achieve a
real DNS for case 1. However the finest grid used to simulate case 1 had only
$512^{3}=1.35\cdot 10^{8}$ grid points. This is an indication that the maximum
grid resolution used in this work is not sufficient.  In addition, as we are
considering two-phase flows with severe jumps in the physical properties,
especially viscosity, properly capturing the interfacial vortical layers
requires that a sufficient number of grid points be located around the
interfaces. The 2D slice shown in Figure 10 clearly indicates that the
vorticity magnitude reaches its maximum near the interfaces, a region where
high shear rates take place, owing to the large viscosity contrast. Therefore
we may suspect that the way the local viscosity is estimated as a function of
the local volume fraction plays a role in the local vorticity magnitude.
Indeed, an arithmetic average is generally performed but there are good
reasons to rather favor an harmonic average [81, 82, 83]. To avoid this
possible influence, we re-computed case 1 with the DyJeAT code while
considering the same viscosity in both fluids. The corresponding results are
discussed in Appendix B. They show that grid convergence of enstrophy is still
not achieved when the viscosity is set to $0.1Pa.s$ (the corresponding
Reynolds number is 130) although the enstrophy excursions are less violent. A
similar study was performed by some of us [84] with the conclusion that sheet
breakup may be responsible for spurious large enstrophy even at reduced Re and
La. The obvious conclusion of these additional computations is that the
differences found among the enstrophy evolutions are not due to effects of an
imprecise numerical modeling of the viscosity jump but to other effects
requiring further study. We also note that in the atomization simulations of
some of us [12] enstrophy as measured statistically is not diverging upon grid
refinement. Hence the results produced by the five codes do not correspond to
a Direct Numerical Simulation: convergence is achieved on quantities dominated
by the large-scale motions, such as the kinetic and potential energies, but it
is not on enstrophy for which the main contribution is from the small-scales.

###  4.2 Case 2: a phase inversion problem involving many break-up events

Fluid | Maximum Enstrophy
---|---
1 |  $103$ $m^{3}.s^{-2}$
2 |  $1049$ $m^{3}.s^{-2}$
Table 6: Values of enstrophy used to make the plots nondimensional in case 2
compared to 5, obtained from DyJeAT on a $2048^{3}$ grid.

In this section, we comment on the simulations of case 2 whose specific
characteristics are provided in Table 2. The corresponding dimensionless data
are given in Table 3. The characteristic time is now
$t_{c,2}=H_{2}/(2U_{g,2})=2.032$ and the dimensionless time is now
$t^{*}=t/t_{c,2}$.

Again, in order to examine grid convergence, the numerical simulations were
run on five grids i.e. $128^{3}$, $256^{3}$, $512^{3}$, $1024^{3}$ and
$2048^{3}$ for DyJeAT. (It has to be noted that 30 million hours of CPU time
have been necessary during two years to perform the simulation on the finest
$2048^{3}$ grid.) This study is displayed on Figure 11. The data of the three
finest grids superpose until approximately $t^{*}=4.5$. After that time, the
comparison of the three finest grids shows that there is no longer a well-
characterized convergence: the predicted kinetic energy varies irregularly
with the grid size, with the difference between the $2048^{3}$ and the
$1024^{3}$ often larger than the difference between the $1024^{3}$ and the
$512^{3}$ grid. This behavior starts at $t^{*}=1.75$ and becomes marked after
$t^{*}=5$. This can only indicate one thing: despite the kinetic energy being
mostly in the large scales, the small scales influence strongly the evolution
of the kinetic energy. These small scales may for example be related to
coalescence events involving the perforation of thin liquid sheets, which will
have a strong influence on the large scales.

Because of the observed convergence behavior of DyJeAT, its $2048^{3}$ kinetic
energy can be taken as a reference, with certainty until $t^{*}=1.75$, and
roughly speaking until $t^{*}=5$. The other numerical simulations were run on
a $256^{3}$ grid for Archer, Jadim, Gerris and Thétis. Similar to the
observations in case 1, the results provided by the five codes are in
approximate agreement as far as the time evolution of the potential and
kinetic energies is concerned, even if larger discrepancies are observed
compared to case 1. These evolutions are displayed in Figures 12 and 13. The
behaviors of the codes are diverse. For example, the Gerris code performs well
for the kinetic energy of fluid 2. However for the kinetic energy of fluid 1
the Gerris code predicts a significantly different oscillatory dynamics as
soon as $t^{*}=1.5$. In fact, for fluid 1, all codes are close to the
reference, except Gerris, until time $t^{*}=2.5$ and then all codes diverge
from the reference. For fluid 2 Gerris and Thetis behave better than the other
codes at times $2<t^{*}<3$.  We have no explanation for this behavior.

An oscillating behavior is still present, even though it it less regular than
in case 1. In all simulations, the kinetic energy decays as $t^{*-2}$ in fluid
2, as already found with case 1. In fluid 1, the exponential Stokes decay law
is clearly more difficult to obtain as large amplitude variations are
observed. Compared to case 1, the turbulent intensity in fluid 2 is larger,
especially in the bottom part of the box, and thus provides a strong forcing
to the large blob of light fluid that stands on top of it. This is why the
oscillatory motion of fluid 1 is more complex than in case 1 and its decay
does not strictly follow the purely viscous Stokes law.

Figure 11: Grid convergence on kinetic energies for case 2 with DyJeAT in log-
linear (top) and log-log (bottom) coordinates.

Figure 12: Potential energies for case 2.

Figure 13: Kinetic energies for case 2 in log-linear (top) and log-log
(bottom) coordinates.

Figure 14: Enstrophy for case 2. Figure 15: Enstrophy for case 2 in log-log
coordinates.

Figure 16: Convergence of enstrophy for case 2.

The time histories of enstrophy are plotted in Figures 14, 15 and 16. Again,
all codes find the peak value nearly at the same time, namely $t\approx 2.5$
in fluid 1 and $t^{*}\approx 4$ in fluid 2, although the code-to-code
differences are significantly larger than in case 1. The delay in the
development of enstrophy in fluid 2 (compared to fluid 1) was not observed in
case 1. It is a direct consequence of the higher Reynolds number: vortical
layers develop quite quickly in the more viscous fluid 1, but a significantly
longer time is required for shear regions to develop in fluid 2, owing to the
prevalence of inertial effects. Again, the enstrophy magnitude increases with
the grid resolution and is code-dependent. Compared to case 1, this tendency
is reinforced by the large number of break-up events which result in a large
population of droplets of fluid 1 that modulates the motion in fluid 2.
Whatever the grid resolution, including $2048^{3}$, convergence is not
achieved and the finer the grid, the larger the enstrophy magnitude. As Figure
16 shows, the difference in magnitude between the peak magnitudes obtained
with two successive resolutions increases as the grid is further refined. This
is an indication that the highest resolution considered here is still far from
that required to achieve a true DNS. However, in the DyJeAT convergence study,
it is observed that the width of the time interval in which the enstrophy
differs is decreasing with refinement of the grid, testifying that grid
convergence will be reached for a finer mesh. This observation is mainly valid
for fluid 1.

Figure 17: Topology of interfaces for case 2 at $t^{*}=4$ (corresponding to
the enstrophy peak) in fluid 2 - Results are from Thétis on a $256^{3}$ grid.

Figure 18: Interfacial area for case 2 using either equation (7) or in the
special case of Gerris the PLIC reconstruction. Equation (7) is noted Int.
Modemi.

A view of the shape of the interfaces at the time the enstrophy is maximum is
provided in Figure 19. Large interfacial scales are observed in the top part
of the box, with at the same time an important population of small dispersed
droplets of fluid 1. Consequently, high-shear regions are more numerous and
difficult to capture than in case 1. The vorticity in fluid 2 is modulated by
the presence of the light droplets of fluid 1 that first follow the motion of
fluid 2 at short times, owing to inertia and viscous effects, and then rise
under buoyancy effects. The finer the grid, the smaller the droplets captured
by the computation. This feature has a direct influence on the intensity and
density of vortical regions in the flow, since the latter are closely related
to the phase distribution, i.e. to the dispersion of fluid 1 in the present
case. Therefore the behavior of the volume-averaged enstrophy depicted above
is a direct consequence of this grid-dependent population of droplets.

In Figure 18 we consider the total interfacial area. For reasons of
simplicity, we use the rough approximation for the interfacial area given in
equation (7). A convergence study for formula (7) is proposed in Appendix C in
the case of a sphere. It is demonstrated that as soon as the grid is refined
the error on the interfacial area is going from $50\%$ to less than $20\%$.

Consider first simulations run on the five grids of DyJeAT. The predicted
interfacial area seems to be almost converged on the finest grid.
Interestingly the convergence is not monotonic, but the interfacial area first
increases for grids $128^{3}$ to $1024^{3}$ then decreases slightly for the
$2048^{3}$ grid. This may indicate that there are two types of error involved,
one coming from the ability of the scheme to capture correctly all the
relevant droplet sizes, and the second from the ability of the code to find,
using (7), the area of a droplet that has been correctly formed, at the
converged size in DyJeAT. The first effect leads to a larger population of
small droplets being produced as the grid is refined from $128^{3}$ to
$1024^{3}$ while the second effect presumably comes from a weaker increase of
the number of droplets but a better capture of their area. One can also see
that in the Gerris case, where the PLIC surface area is available, the
interfacial area is markedly smaller with PLIC than with (7).

The physical mechanism behind the observed interfacial area behavior is easy
to comprehend in its broad lines. The initial fast upward motion of fluid 1
atomizes the interface and results in the growth of a population of small
droplets, that yields an increase of the interfacial area. However,
mysteriously, it can be noticed that the peak value of the interfacial area is
not only grid-dependent but also code-dependent. Pure PLIC-VOF codes such as
Gerris conserve mass much better and the breakup of thin sheets results in the
formation of a large number of droplets, while level-set codes such as DyJeat
result in the “evaporation” of small objects. Thus the relatively larger
number of droplets at late times in Jadim and Gerris, both VOF codes, and the
smaller interfacial area in DyJeAT, a Level-Set code, could be attributed to
this fact and not to the use of (7) to compute the interfacial area. However
this does not explain why Thetis, also a PLIC-VOF code, has such a small
interface area. Other aspects of the codes may explain the properties of
droplet atomization and coalescence, such as the implementation of surface
tension forces, which varies widely between codes.

This issue needs to be addressed in future investigations because the
interfacial area is one of the key parameters in problems involving heat or
mass transfer.

##  5 Droplet size distributions

We have performed a study of droplet size distributions using DyJeAT in case 2
at the time of the enstrophy peak ($t^{*}=4$). We define the Probability
Distribution Function (PDF) as number frequency $f$ defined as follows. The
droplet sizes are measurered in the simulations and an equivalent radius $r$
is obtained for each drop. The equivalent radius may belong in any of a number
of adjacent bins $B_{i}$ defined as the interval
$B_{i}=(r_{i}-\Delta_{i}/2,r_{i}+\Delta_{i}/2)$ where $r_{i}$ and $\Delta_{i}$
are the bin center and width, and define the bin structure. In our first
analysis, the bins are regularly spaced. One defines $\Delta_{i}=1/256$,
$r_{i}=(i-1/2)/256$ and $B_{i}=((i-1)/256,i/256)$. Then the number of droplets
whose equivalent radius is in $B_{i}$ is noted $N_{i}$. $N_{i}$ is
approximately proportional to an asymptotic probability $f$ such that
$f(r_{i})=N_{i}/\Delta_{i}$ however the approximation does not hold very well
for the first bin because $\Delta_{1}/r_{1}$ is far from small. On Figure 19
we plot $N_{i}$ versus $r_{i}$.

Figure 19: Droplet radius frequency $N_{i}$ defined in the text obtained with
DyJeAT at $t^{*}=4$. The droplet radii are in meters and the box size is
$H=1$m.

The observation of these distributions leads to several remarks. First, the
$256^{3}$ grid result does not match the others, but instead droplet
populations obtained with that grid are systematically less numerous than for
the other grids, indicating that the total mass of droplets obtained with that
grid is smaller. This can be easily explained by the fact that poorly resolved
droplets “evaporate” in a level-set method such as DyJeAT. The uncertainty
surrounding the $256^{3}$ grid results also leads us not to discuss the
results of the other codes for the PDF, which are only available at this
resolution. Second, the $512^{3}$, $1024^{3}$ and $2048^{3}$ grid results
superpose in a region of variable size, going from a minimum radius $r_{m}(N)$
to a maximum radius $r_{M}(N)$, with both $r_{m}$ and $r_{M}$ function of the
grid resolution $N$. To be specific, the minimum radius of agreement decreases
with increasing $N$ until approximately $r_{m}(2048)\simeq 0.005$ which
indicates a range of convergence extending to very small droplet radii for the
$2048^{3}$ grid. The upper boundary $r_{M}(N)$ of the range of convergence is
harder to determine, because the large radii are affected by statistical
noise. However, it is seen that few additional very large droplets ($r>0.07$)
are seen as the resolution increases from $256^{3}$ to $2048^{3}$. Moreover
grids $1024^{3}$ and $2048^{3}$ are in agreement in a wide range with
$r_{M}(2048)\simeq 0.05$.  Finally we note that the behavior at very small
droplet sizes is different for the $2048^{3}$ grid and other grids: while the
other grids have a cutoff (a fall in droplet counts $N_{1}$ and/or $N_{2}$) at
small sizes, the $2048^{3}$ grid does not. This is due to the fact that in the
$2048^{3}$ case the fall in droplet count can only be seen if instead of the
peculiar bin $B_{1}$ defined above, one subdivides the interval $(0,1/256)$
into several smaller bins. We have actually done that and did observe the
cutoff even for the $2048^{3}$ case.

Figure 20: Droplet radius frequencies (probability distribution function) in
logarithmic coordinates, with $rf(r)$ versus $r$. A log-normal distribution
and an exponential one are shown to guide the eye.

Indeed, similar effects of the grid refinement, including the cutoff or drop
in frequency, have been reported for atomization of liquid jets, especially
for level set methods [85] and less markedly for VOF methods [19]. This
indicates that $512^{3}$ grids are a minimum requirement to have an
approximately accurate droplet count at any scale, and that $256^{3}$
simulations are poorly resolved in the whole range of scales. This explains
the poor results on interfacial area with the $256^{3}$ simulations of all the
other codes. On the other hand, the absence of convergence of the enstrophy
does not imply a non convergence of the droplet sizes in the intermediate
range $r_{m}(N)<r<r_{M}(N)$. Another interesting fact is that the large scales
are influenced by how the small scales are computed: more large droplets are
seen with finer grids, probably because more grid resolution implies less
breakup of thin layers or filaments, thus reducing the rate at which large
structures break into small ones.  In order to better understand the PDF, we
plot the computed frequencies for the two most refined grids as a graph of
$\ln(r_{i}N_{i})$ versus $\ln r_{i}$ in Figure 20. The first bin $B_{1}$ has
been excluded due to its particular character, and only the most refined grids
$1024^{3}$ and $2048^{3}$ are shown. Using these coordinates has the advantage
that the log-normal PDF appears as a parabola. The log-normal PDF is defined
by

$f(r)=\frac{A}{r}\exp\left[-\frac{(\ln
r-\ln\hat{\mu})^{2}}{2\hat{\sigma}^{2}}\right]$ (12)

where $A$ is normalisation constant, $\ln\hat{\mu}$ is the logarithmic average
and $\hat{\sigma}$ the logarithmic standard deviation. A parabolic curve
corresponding to such a log-normal distribution is shown on Figure 20 to aid
the interpretation. This parabola is built with $\hat{\mu}=10^{-3}m$,
$\hat{\sigma}=1.7$ and $A=3$. While the value of $\hat{\mu}$ is at the least
very uncertain, it is clear that extrapolating the trend of the distribution
leads to expect a larger value of $rf(r)$ at smaller $r$. These larger values
of $rf(r)$ could be revealed by simulations with yet more refined grids than
the current ones. This strongly reinforces the expectation that many smaller
droplets would be seen in a true DNS. However these droplets would contribute
little to the interface area. Indeed the interface area scales as

$S=\int_{0}^{\infty}4\pi r^{2}f(r)\,dr$ (13)

which for a log-normal $f(r)$ or for a simpler $f(r)\sim 1/r$ converges at the
$r=0$ bound. Note that we took the assumption of spherical droplets to
estimate $S$. The latter assumption is likely to be a good approximation for
the smallest droplets. Thus the absence of the smallest droplets in the
integral does not affect much the interface area $S$.

A simpler distribution of the form $f(r)\sim 1/r$ would be a horizontal line
in the variables of Figure 20. That approximation would not be much worse than
the log-normal one in the intermediate range of droplet sizes. As a third
alternative, an exponential distribution $f(r)=A\exp(-r/r_{m})$ is also
plotted on Figure 20 with $r_{m}=0,02$ m. The exponential distribution is a
better fit for the largest values of $r$. With this distribution, the
interfacial area integral converges for $\Delta x\ll r_{m}=0,02\,$m and thus
for $256^{3}$ grids. The fact that the interfacial area is seen not to
converge in simulations signifies that the number of droplets is much below
its theoretical level, indicating a kind of numerical “evaporation” of the
droplets.

Finally we note that none of the distributions yields a very good fit. It is
possible that the PDF is in fact bimodal, with different “modes” being
superposed at various scales. This could be clarified by computing on even
more refined grids to capture the small sizes, or using an ensemble of initial
conditions as discussed below.

##  6 Summary and future work

A computational benchmark based on phase inversion of two immiscible fluids of
different densities confined in a closed cubic box has been set up. The
purpose of this benchmark is to check the capabilities and limitations of
current codes and grid resolutions with the ultimate goal of producing
numerical data of a quality similar to those provided by a DNS of single-phase
turbulence. Several numerical approaches based on the one-fluid model and
routinely used to compute two-phase incompressible flows have been used and
compared. They include the Volume-Of-Fluid and Level-Set approaches for
interface tracking, the Ghost-Fluid technique for the capture of interfacial
jumps of mass and momentum, as well as several variants of Navier-Stokes
solvers such as the time-splitting and augmented-Lagrangian algorithms. All
these approaches and methods are a priori suitable for simulating the various
configurations encountered in phase-inversion problems on an arbitrary grid.

The first problem we considered, case 1, is characterized by an integral scale
Reynolds number Re = $1.37\,10^{4}$. Computational results have been analyzed
by considering the time evolution of several volume-averaged indicators,
namely the potential and kinetic energy in each fluid, the relative volume of
light fluid in the top part of the box and the enstrophy in each fluid. The
results reveal that all codes provide close evolutions for the first three
quantities, indicating that reliable predictions are obtained when quantities
essentially determined by large-scale motions are considered. The long-term
decay rates of the kinetic energies have been found to be in good agreement
with the viscous Stokes law (in fluid 1) and the turbulent decay law (in fluid
2) discussed in Appendix A, respectively. In addition, all codes agree on the
sloshing frequency of the light fluid. The situation has been found to be much
more problematic when enstrophy is considered: all codes provide markedly
different values of the enstrophy maximum (although they agree on the time at
which this maximum occurs) and none of them converges with respect to grid
refinement.

A second phase inversion test case, involving a higher Reynolds number Re =
$4.33\,10^{5}$was subsequently considered. The same general conclusions apply
to this configuration. All codes provide converging predictions as long as
only kinetic and potential energies are considered, whereas large
discrepancies are observed on enstrophy and interfacial area. Here again, the
DNS conditions are not satisfied and much finer grids should be considered to
expect convergence on quantities such as enstrophy and interfacial area which
are governed by small-scale processes, especially multiple break-up events. It
has to be noticed that despite the discrepancies observed on enstrophy,
approximate convergence is obtained for PDFs of droplet sizes as soon as a
$1024^{3}$ grid is used. Clearly, case 2 is an implicit LES rather than a DNS.
As a conclusion, this benchmark illustrates the fact that implicit LES
simulations on fine grids allow to reach convergence on first order moments
such as kinetic or potential energies while they fail to converge as soon as
small scale statistics are investigated. Explicit LES models may be
implemented in this case [34, 35, 24, 18, 86, 87] or finer meshes should be
utilized for true DNS.

In order to disentangle physically meaningful phenomena and numerical
artifacts possibly induced in the interfacial regions by the averaging rule
used to define the local viscosity, complementary problems have been
considered in Appendix B. They make use of the same parameters as case 1,
except for viscosity $\mu$ which is the same in both fluids and equal to
$0.1$, $0.01$ and $0.001Pa.s$. For the two largest viscosities ($0.1$ and
$0.01Pa.s$), convergence is still not obtained for enstrophy in both fluids,
but less violent excursions are observed.

Finally, a natural and important extension of this study should address the
reproducibility of the benchmark results. In the approach reported in the
current paper, we do not select or force any initial interface perturbation by
specifying perfectly flat interfaces located on the faces of a cube. No
initial wavelength or length scale other than the size of the cube is present.
Constitently with this flat character of the interface the initial velocity is
zero. In an extension of this study, we have begun investigating a statistical
ensemble of simulations. The ensemble is constructed by adding an interfacial
or velocity perturbation at $t=0$. We then perform repeatedly the same
simulation with a different random perturbation each time. A statistical
analysis of the problem is then performed in order to extract an average or
typical flow field and interface position, and fluctuations about this
average. This approach is in progress. Clearly, characterizing how simulations
that are part of an ensemble with noisy or perturbed initial conditions
diverge from each other is an important issue that needs to be investigated
and constitutes an important perspective of extension for the present study.

Future work of all partners will first be devoted to reconsidering cases 1 and
2 on finer grids. Multiscale Eulerian VOF or Level Set approaches coupled to a
Lagrangian description of small droplets will also be introduced by some of
them with the goal of reaching convergence on all physical quantities on
‘reasonable’ grids that do not capture the smallest droplets. Last, adaptative
mesh refinement (AMR) techniques will be considered to concentrate numerical
efforts in flow regions where the vorticity magnitude is expected to reach its
peak.

##  7 Acknowledgements

This work is part of the MODEMI ANR project (Modélisation et Simulation Multi-
échelle des Interfaces, ANR-11-MONU-0011) devoted to the multiscale modeling
of two-phase flows. We thus acknowledge the support of the Agence Nationale de
la Recherche (ANR) for the I2M, DALEMBERT, CORIA and IMFT laboratories
concerning the Multiscale Modeling of Interfaces. The authors wish to thank
the Midi-Pyrénées and Aquitaine Regional Councils for the financial support
dedicated to a PhD thesis at ONERA and IMFT and a 256-processor cluster
investment, located in the TREFLE team of the I2M laboratory. We are grateful
for access to the computational facilities of the French CINES (National
computing center for higher education) and CCRT (National computing center of
CEA) under project number x20142b6115, x20152b6115, x20162b6115 and
A0012b06115. This work was also partly granted access to the HPC resources of
CALMIP under the allocation 2013-P0633. The CORIA team would like to express
its gratitude to the CRIHAN (Centre de Ressources Informatiques de HAute-
Normandie, www.crihan.fr) computational centre used for providing CPU
resources.

##  Appendix A Appendix A: Scaling laws for the kinetic energy decay

After the acceleration stage induced by buoyancy forces, during which most of
the lighter fluid goes to the top part of the cavity, the phase inversion
problem is characterized in a second time by a wavy behavior which makes it
look quite similar to a sloshing flow progressively damped by shearing and
viscous effects. The classical analysis of Stokes regarding the viscous
damping of gravity waves (see e.g. [88], pp. 623-624), makes it possible to
predict the time evolution of the kinetic energy in a weakly viscous flow
driven by a surface wave. First, the time evolution of the mechanical energy,
which is the sum of the potential and kinetic energies, is known to result
from the internal viscous dissipation, so that

$\frac{dE_{mech}}{dt}=-\frac{1}{2}\mu\int_{\Omega}\left(\nabla{\bf
u}+\nabla^{t}{\bf u}\right)^{2}dV$ (14)

If we assume the shear layer at the free surface to have negligible effect
owing to the moderate velocity gradients expected in this region, the flow can
be considered as irrotational. Then the velocity potential, $\phi$, and
mechanical energy, $E_{mech}$, read

$\displaystyle\phi=\phi_{0}e^{kz}cos(kx-\omega t)\,,$ (15)
$\displaystyle\frac{d\overline{E_{mech}}}{dt}=-8\mu
k^{4}\int_{\Omega}\overline{\phi^{2}}dV\,,$ (16)

where $k$ and $\omega$ are the wave number and wave frequency, respectively,
and the overbar denotes the time average value. In the framework of linear
wave theory, the potential and kinetic energies are equal. However the wave
amplitude is not small in the present phase inversion problem. Nevertheless,
as soon as most light fluid lies in the top part of the cavity, the potential
energy stays almost constant, making the time variation of $E_{mech}$ mainly
governed by that of the kinetic energy. Therefore, still in the linear
approximation, we can approximately write

$\overline{E_{mech}}=\frac{1}{2}\rho\int_{\Omega}\overline{{\bf
u}^{2}}dV\approx\frac{1}{2}\rho\int_{\Omega}{\overline{\nabla\phi\cdot\nabla\phi}}dV=\rho
k^{2}\int_{\Omega}\overline{\phi^{2}}dV$ (17)

Now, two markedly different flow situations can be met in the phase inversion
problem,namely a “gentle” configuration in which break-up and coalescence
events of the light fluid is scarcely occur and a “violent” configuration in
which such events are much more numerous. Combining (16) and (17), the time
evolution of the kinetic energy in the first regime is found to obey

$\overline{E_{mech}}=\frac{1}{2}\rho\int_{\Omega}\overline{{\bf
u}^{2}}dV=Ke^{-8\nu k^{2}t}=Ke^{-8\nu\frac{\omega^{4}}{g^{2}}t}\,,$ (18)

where $K$ is a constant, $\nu=\mu/\rho$, and use has been made of the
dispersion relation $\omega^{2}=gk$.
Let us now consider the “violent” configuration. In this case, averaging
throughout the whole volume $\Omega$, (15)-(16) simply yields

$\displaystyle\widetilde{{\bf
u}^{2}}=\widetilde{\nabla\phi\cdot\nabla\phi}=k^{2}\widetilde{\phi^{2}}\,,$
(19) $\displaystyle\frac{d\overline{E_{mech}}}{dt}=-8\mu k^{2}\widetilde{{\bf
u}^{2}}\,,$ (20)

where $\widetilde{\cdot}$ denotes the mean value over $\Omega$. In other
words, one has

$\frac{d\widetilde{{\bf u}^{2}}}{dt}=-16\nu k^{2}\widetilde{{\bf u}^{2}}$ (21)

In this configuration, the flow is expected to be turbulent and the molecular
viscosity is not relevant any more to estimate the flow damping rate. As a
crude surrogate, we can introduce an effective turbulent viscosity $\nu_{t}$
scaling as $l\widetilde{\bf u}$, where $l$ is a characteristic length of the
large-scale flow, which can be taken as the box size. Replacing $\nu$ by
$\nu_{t}$ in (21), we then obtain

$\frac{d\widetilde{{\bf u}^{2}}}{dt}=-16lk^{2}\widetilde{{\bf u}^{3}}$ (22)

Assuming the volume-averaged velocity to follow a power law, i.e.
$\widetilde{{\bf u}}$ to be of the form $\widetilde{{\bf u}_{0}}t^{-n}$,
immediately yields $n=1$, from which we infer

$\overline{E_{mech}}\sim t^{-2}$ (23)

This prediction is to be compared with the exponential decay predicted by (18)
in the “gentle” regime. The result (23) is reminiscent of the approximate
decay law of decaying homogeneous isotropic turbulence (e.g. [89]). This is
expected since, according to Taylor’s estimate, the dissipation rate
$.pngilon$ is known to scale as $u_{0}^{3}/l_{0}$, where $u_{0}$ and $l_{0}$
stand for the large-scale velocity and length scales, respectively. Therefore,
when $l_{0}$ is constant, (23) is immediately recovered and the dissipation
rate is predicted to decay as $t^{-3}$.

##  Appendix B Appendix B: influence of viscosity on the flow dynamics

In order to better understand the origin of code- and grid-dependences of the
volume-averaged enstrophy in cases 1 and 2, we perform an extra series of
computations in which the viscosity is identical in the two fluids. In this
way, any possible influence of the averaging procedure selected to compute the
local viscosity as a function of the volume fraction and of the numerical
treatment of the viscosity jump in the interfacial grid cells is removed. We
make use of the physical parameters of case 1, except for viscosity which is
set to $0.1Pa.s$ (case 1a), $0.01Pa.s$ (case 1b) and $0.001Pa.s$ (case 1c),
respectively. The simulations are only carried out with the DyJeAT code, as it
allows us to use the largest grids at a reasonable cost. The convergence study
is performed on three different grids, namely $128^{3}$, $256^{3}$ and
$512^{3}$ in case 1a, whereas four grids ranging from $128^{3}$ to $1024^{3}$
are considered in cases 1b and 1c.

Figure 21: Grid convergence of enstrophy for case 1a in fluid 1 (left) and
fluid 2 (right) - linear-linear (top) and log-linear (bottom) coordinates.

Figure 22: Grid convergence of enstrophy for case 1b in fluid 1 (left) and
fluid 2 (right) - linear-linear (top) and log-log (bottom) coordinates.

Figure 23: Grid convergence of enstrophy for case 1c in fluid 1 (left) and
fluid 2 (right) - linear-linear (top) and log-log (bottom) coordinates.

The evolution of the volume-averaged enstrophy is reported in Figures 21, 22
and 23 for cases 1a, 1b and 1c, respectively. In case 1a, the Reynolds number
is $130$ and the enstrophy while still not converging as the grid is further
refined, has less violent excursions. In this configuration, the enstrophy
decay at large time follows an exponential law, a behavior typical of the
viscous decay of gravity waves (see appendix A). In case 1b, grid convergence
is also not achieved in both fluids, the Reynolds number in fluid 2 being
$1300$. At large time, the enstrophy obeys a $t^{*-3}$ decay law typical of
turbulent conditions (see appendix A). In case 1c, enstrophy convergence is
not achieved in either fluid. The Reynolds number is $1.37\,10^{4}$and it is
observed that, even though the results concerning energy and relative volume
in the upper part of the box converge, as shown with case 1, the grid is not
thin enough to capture the small-scale eddies. Again, the decay law at large
times corresponds of course to turbulent conditions.

##  Appendix C Appendix C: convergence of the interfacial area

A simple formula is proposed in the present work to estimate the interfacial
area. It is based on the detection of the cells cut by the interface by
considering the gradient of the volume of fluid. It is reported in equation
(7). This estimate can be considered to be rough compared to a more accurate
calculation provided by the sum of the surface of the planes generated by the
VOF-PLIC method. However, the interest of (7) is that it can be utilized
whatever the interface tracking algorithm (VOF or Level Set).

Figure 24: Convergence of formula (7) for interface area calculation - sphere
surface on a $128^{3}$ grid (left) and grid convergence of the interfacial
area (right) compared to the exact sphere surface.

In order to evaluate the accuracy of formula (7), it is proposed to estimate
the interfacial area of a sphere of radius $0.25m$ placed in a square domain
of side size $1m$. The sphere surface obtained with a VOF function on a
$128^{3}$ grid is presented in Fig. 24 together with the convergence of the
percentage error on the interfacial area with respect to grid refinement. It
has to be noticed that in the present case, the exact interfacial area is
$\pi/4$. The percentage error decreases from $50\%$ to $18\%$ for a number of
grid cells being $8$ to $128$ in each Cartesian space direction. The coarser
grid is representative of an under resolved drop containing $2$ grid cells
over a radius whereas the $128^{3}$ grid is representative of a well resolved
drop with $32$ grid cells over a radius. These are typical characteristics
that are encountered in the phase inversion simulations. It can be estimated
that the error brought by (7) is belonging to the range $20\%$ to $50\%$ which
is not so bad for a rough approximation. A better estimate is certainly
provided by the sum of the surface of the planes coming from the VOF-PLIC
reconstruction. An illustration is given in Fig. (18). The gap between (7)
(noted Int. Modemi) and VOF PLIC calculation (noted Int. PLIC) on a $128^{3}$
grid for the Gerris code is $30\%$, which is in agreement with the present
convergence study.

## References

  * [1] W.J. Rider and D.B. Kothe.  Reconstructing volume tracking.  J. Comput. Phys., 141:112–152, 1998.
  * [2] M. Rudman.  A volume-tracking method for incompressible multifluid flows with large density variations.  Int. J. Numer. Meth. Fluids, 28:357–378, 1998.
  * [3] Various Authors.  Test-cases for interface tracking methods.  http://test.interface.free.fr/, 2003.
  * [4] A. Benkenida and J. Magnaudet.  A method for the simulation of two-phase flows without interface reconstruction.  C. R. Acad. Sci. Paris Série IIb, 328:25–32, 2000.
  * [5] S. Popinet.  An accurate adaptive solver for surface-tension-driven interfacial flows.  J. Comput. Phys., 228:5838–5866, 2009.
  * [6] A. Prosperetti.  Motion of two superposed viscous fluids.  Phys. Fluids, 24:1217, 1981.
  * [7] N. Ashgriz and J.Y. Poo.  Coalescence and separation in binary collisions of liquid drops.  J. Fluid Mech., 221:183–204, 1990.
  * [8] S. Tanguy and A. Berlemont.  Application of a level set method for simulation of droplet collisions.  Int. J. Multiphase Flow, 31:1015–1035, 2005.
  * [9] J.C. Martin and W.J. Moyce.  An experimental study of the collapse of liquid columns on a rigid horizontal plane.  Philos. Trans. R. Soc. A-Math. Phys. Eng. Sci., 244:312–324, 1952\.
  * [10] J. Delteil, S. Vincent, A. Erriguible, and P. Subra-Paternault.  Numerical investigations in Rayleigh breakup of round liquid jets with VOF methods.  Comput. Fluids, 50:10–23, 2011.
  * [11] J. C. Lasheras and E. J. Hopfinger.  Liquid jet instability and atomization in a coaxial gas stream.  Annu. Rev. Fluid Mech., 32:275–308, 2000.
  * [12] Y. Ling, D. Fuster, G. Tryggvason, and S. Zaleski.  A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability.  J. Fluid Mech., 859:268–307, 2019.
  * [13] J. Eggers and E. Villermaux.  Physics of liquid jets.  Rep. Prog. Phys., 71:036601, 2008.
  * [14] L. Gilman and E. Baglietto.  A self-consistent, physics-based boiling heat transfer modeling framework for use in computational fluid dynamics.  Int. J. Multiphase Flow, 95:35–53, 2017.
  * [15] J. Lu and G. Tryggvason.  Direct numerical simulations of multifluid flows in a vertical channel undergoing topology changes.  Phys. Rev. Fluids, 3(8):084401, 2018.
  * [16] C. Josserand and S.T. Thoroddsen.  Drop impact on a solid surface.  Annu. Rev. Fluid Mech., 48:365–391, 2016.
  * [17] L. Deike, W.K. Melville, and S. Popinet.  Air entrainment and bubble statistics in breaking waves.  J. Fluid Mech., 801:91–129, July 2016.
  * [18] S. Vincent, D. Lacanette, J. Larocque, A. Toutant, P. Lubin, and P. Sagaut.  Direct numerical simulation of phase separation and a priori two-phase LES filtering.  Comput. Fluids, 37:898–906, 2008.
  * [19] Y. Ling, D. Fuster, S. Zaleski, and G. Tryggvason.  Spray formation in a quasiplanar gas–liquid mixing layer at moderate density ratios: A numerical closeup.  Phys. Rev. Fluids, 2(1):014005, January 2017.
  * [20] J.-P. Matas, A. Delon, and A. Cartellier.  Shear instability of an axisymmetric air–water coaxial jet.  J. Fluid Mech., 843:575–600, March 2018.
  * [21] M. Lance and J. Bataille.  Turbulence in the liquid phase of a uniform bubbly air-water flow.  J. Fluid Mech., 222:95–118, 1991.
  * [22] H. Chanson.  Air-water flow measurements with intrusive, phase-detection probes: Can we improve their interpretation?  J. Hydraulic Eng., pages 252–255, 2002.
  * [23] K.M.T. Kleefsman, G. Fekken, A.E.P. Veldman, B. Iwanowski, and B.A. Buchner.  Volume-of-fluid based simulation method for wave impact problems.  J. Comput. Phys., 206:363–393, 2005.
  * [24] J. Larocque, S. Vincent, P. Lubin, D. Lacanette, and J.-P. Caltagirone.  Parametric study of LES subgrid terms in a turbulent phase separation flow.  Int. J. Heat Fluid Flow, 31:536–544, 2010.
  * [25] A.W. Cook, W. Cabot, and P.L. Miller.  The mixing transition in Rayleigh-Taylor instability.  J. Fluid Mech., 511:333–362, 2004.
  * [26] P. Ramaprabhu, G. Dimonte, and M.J. Andrews.  A numerical study of the influence of initial perturbations on the turbulent Rayleigh-Taylor instability.  J. Fluid Mech., 536:285–319, 2005.
  * [27] G. Boffetta, F. De Lillo, A. Mazzino, and L. Vozella.  The ultimate state of thermal convection in Rayleigh-Taylor turbulence.  Physica D, 241:137–140, 2012.
  * [28] S. I. Abarzhi, A. Gorobets, and K. R. Sreenivasan.  Rayleigh-Taylor turbulent mixing of immiscible, miscible and stratified fluids.  Phys. Fluid, 17:081705:1–4, 2005.
  * [29] G. Dimonte, D. L. Youngs, A. Dimits, S. Weber, M. Marinak, S. Wunsch, C. Garasi, A. Robinson, M. J. Andrews, P. Ramaprabhu, A. C. Calder, B. Fryxell, J. Biello, L. Dursi, P. MacNeice, K. Olson, P. Ricker, R. Rosner, F. Timmes, H. Tufo, Y.-N. Young, and M. Zingale.  Rayleigh-Taylor turbulent mixing of immiscible, miscible and stratified fluids.  Phys. Fluid, 16:1668–1693, 2004.
  * [30] G. Tryggvason.  Numerical simulations of the Rayleigh-Taylor instability.  J. Comput. Phys, 75:253–282, 1988.
  * [31] W. Aniszewski, T. Ménard, and M. Marek.  Volume of fluid (VOF) type advection methods in two-phase flow: A comparative study.  Comput. Fluids, 97:52 – 73, 2014.
  * [32] D. Jacqmin.  Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling.  J. Comput. Phys, 155:96–127, 1999.
  * [33] S. Vincent, J.-P. Caltagirone, and D. Jamet.  Test case n°15: phase inversion in a closed box (pn, pe).  Multiphase Sci. Technol., 6:101–104, 2004.
  * [34] W. Aniszewski, A. Boguslawski, M. Marek, and A. Tyliszczak.  A new approach to sub-grid surface tension for LES of two-phase flows.  J. Comput. Phys., 231:7368–7397, 2012.
  * [35] E. Labourasse, D. Lacanette, A. Toutant, P. Lubin, S. Vincent, O. Lebaigue, J.-P. Caltagirone, and P. Sagaut.  Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests.  Int. J. Multiphase Flow, 33:1–39, 2007.
  * [36] R. Lebas, T. Ménard, P.A. Beau, A. Berlemont, and F.X. Demoulin.  Numerical simulation of primary break-up and atomization: DNS and modelling study.  Int. J. Multiphase Flow, 35:247–260, 2009.
  * [37] G. Tomar, D. Fuster, S. Zaleski, and S. Popinet.  Multiscale simulations of primary atomization.  Comput. Fluids, 39:1864–1874, 2010.
  * [38] I. Calmet and J. Magnaudet.  High-schmidt number mass transfer through turbulent gas-liquid interfaces.  Int. J. Heat Fluid Flow, 19:522–532, 1998.
  * [39] S. Banerjee, D. Lakehal, and M. Fulgosi.  Surface divergence models for scalar exchange between turbulent streams.  Int. J. Multiphase Flow, 30:963–977, 2004.
  * [40] J.M. Delhaye.  Jump conditions and entropy sources in two-phase systems. Local instant formulation.  Int. J. Multiphase Flow, 1:395–409, 1974.
  * [41] R. Scardovelli and S. Zaleski.  Direct numerical simulation of free-surface and interfacial flow.  Annu. Rev. Fluid Mech., 31:567–603, 1999.
  * [42] J.U. Brackbill, D.B. Kothe, and C. Zemach.  A continuum method for modeling surface tension.  J. Comput. Phys., 100:335–354, 1992.
  * [43] S. Osher and R. Fedkiw.  Level set methods: An overview and some recent results.  J. Comput. Phys., 169:463–502, 2001.
  * [44] S. Patankar.  Numerical heat transfer and fluid flow.  CRC Press, 1980.
  * [45] A. Poux, S. Glockner, and M. Azaïez.  Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods.  J. Comput. Phys., 230:4011–4027, 2011.
  * [46] S. Vincent, A. Sarthou, J.-P. Caltagirone, F. Sonilhac, P. Février, C. Mignot, and G. Pianet.  Augmented Lagrangian and penalty methods for the simulation of two-phase flows interacting with moving solids. Application to hydroplaning flows interacting with real tire tread patterns.  J. Comput. Phys., 230:956–983, 2011.
  * [47] F. Bertrand, P.A. Tanguy, and F. Thibault.  A three-dimensional fictitious domain method for incompressible fluid flow problems.  Int. J. Numer. Meth. Fluids, 25:719–736, 1997.
  * [48] M. Fortin and R. Glowinski.  Méthodes de lagrangien augmenté. Application à la résolution numérique de problèmes aux limites.  Dunod, Paris, 1982.
  * [49] O. Schenk and K. G$\ddot{a}$rtner.  Solving unsymmetric sparse systems of linear equations with PARDISO.  Fut. Gener. Comput. Syst., 20:475–487, 2004.
  * [50] RH. A. Van Der Vorst.  BiCGSTAB: a fast and smoothly converging varian of BiCG for the solution of non-symmetric linear systems.  J. Sci. Comput., 13:631–644, 1992.
  * [51] D.L. Youngs.  Time-dependent multimaterial flow with large fluid distortion.  K.W. Morton and M.J. Baines eds., Numerical Methods for Fluid Dynamics, Academic, New-York, 1982.
  * [52] S. Vincent, G. Balmigère, J.-P. Caltagirone, and E. Meillot.  Eulerian-Lagrangian multiscale methods for solving scalar equations \- Application to incompressible two-phase flows.  J. Comput. Phys., 229:73–106, 2010.
  * [53] T.N. Randrianarivelo, G. Pianet, S. Vincent, and J.-P. Caltagirone.  Numerical modelling of the solid particle motion using a new penalty method.  Int. J. Numer. Meth. Fluids, 47:1245–1251, 2005.
  * [54] S. Vincent, T. Randrianarivelo, G. Pianet, and J.P. Caltagirone.  Local penalty methods for flows interacting with moving solids at high Reynolds numbers.  Comput. Fluids, 36:902–913, 2007.
  * [55] E. Aulisa, S. Manservisi, R. Scardovelli, and S. Zaleski.  Interface reconstruction with least-squares fit and split advection in three-dimensional Cartesian geometry.  J. Comput. Phys., 225:2301–2319, 2007.
  * [56] J. Li.  Calcul d′interface affine par morceaux (piecewise linear interface calculation).  C. R. Acad. Sci. Paris Série IIb, 320:391–396, 1995.
  * [57] R. Scardovelli and S. Zaleski.  Interface reconstruction with least-square fit and split Lagrangian-Eulerian advection.  Int J. Numer. Meth. Fluids, 41:251–274, 2003.
  * [58] E. Aulisa, S. Manservisi, R. Scardovelli, and S. Zaleski.  A geometrical area-preserving volume of fluid advection method.  J. Comput. Phys., 192:355–364, 2003.
  * [59] G. Tryggvason, R. Scardovelli, and S. Zaleski.  Direct numerical simulations of gas-liquid multiphase flows.  Cambridge University Press, 2011.
  * [60] M.M. Francois, S.J. Cummins, E.D. Dendy, D.B. Kothe, J.M. Sicilian, and M.W. Williams.  A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework.  J. Comput. Phys., 213:141–173, 2006.
  * [61] Y. Renardy and M. Renardy.  PROST: A parabolic reconstruction of surface tension for the volume-of-fluid method.  J. Comput. Phys., 183:400–421, 2002.
  * [62] S. Popinet.  Numerical models of surface tension.  Annu. Rev. Fluid Mech., 50:49–75, 2018.
  * [63] D. Fuster, J.-Ph. Matas, S. Marty, S. Popinet, J. Hoepffner, A. Cartellier, and S. Zaleski.  Instability regimes in the primary breakup region of planar coflowing sheets.  J. Fluid Mech., 736:150–176, 2013.
  * [64] D. Fuster, A. Bagué, T. Boeck, L. Le Moyne, A. Leboissetier, S. Popinet, P. Ray, R. Scardovelli, and S. Zaleski.  Simulation of primary atomization with an octree adaptive mesh refinement and VOF method.  Int. J. Multiphase Flow, 35:550–565, 2009.
  * [65] D. Fuster, G. Agbaglah, C. Josserand, S. Popinet, and S. Zaleski.  Numerical simulation of droplets, bubbles and waves: state of the art.  Fluid Dyn. Res., 41:065001, 2009.
  * [66] T. Bonometti and J. Magnaudet.  An interface-capturing method for incompressible two-phase flows. validation and application to bubble dynamics.  Int. J. Multiphase Flow, 33:109–133, 2007.
  * [67] B. Duret, J. Reveillon, T. Menard, and F.X. Demoulin.  Improving primary atomization modeling through DNS of two-phase flows.  Int. J. Multiphase Flow, 55:130 – 137, 2013.
  * [68] Chi-Wang Shu.  Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws.  NASA (as NASA/CR-97-206253 report), 1997.
  * [69] M. Sussman, E. Fatemi, P. Smereka, and S. Osher.  Improved level set method for incompressible two-phase flows.  Comput. Fluids, 27:663–680, 1998.
  * [70] M. Sussman and E.G. Puckett.  A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows.  J. Comput. Phys., 162:301–337, 2000.
  * [71] M. Kang, R.P. Fedkiw, and X.D. Liu.  A boundary condition capturing method for multiphase incompressible flow.  J. Sci. Comput., 15:323–360, 2000.
  * [72] X.D. Liu, R.P. Fedkiw, and M. Kang.  A boundary condition capturing method for poisson equation on irregular domains.  J. Comput. Phys., 260:151–178, 2000.
  * [73] J. D. Towers.  Finite difference methods for approximating heaviside functions.  J. Comput. Phys., 228:3478–3489, 2009.
  * [74] A.J. Chorin.  Numerical solution of the Navier-Stokes equations.  Math. Comput., 22:745–762, 1968.
  * [75] R. Temam.  Sur l’approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires.  Arch. Ration. Mech. Anal., 33:377–385, 1969.
  * [76] R. Fedkiw, T. Aslam, B. Merriman, and S. Osher.  A non-oscillatory Eulerian approach to interfaces in multimaterials flows (the ghost fluid method).  J. Comput. Phys., 152:457–492, 1999.
  * [77] D. Zuzio and J.-L. Estivalezes.  An efficient block parallel AMR method for two phase flow simulations.  Comput. Fluids, 44:339–357, 2011.
  * [78] D. Zuzio, J.-L. Estivalèzes, and B. DiPierro.  An improved multiscale Eulerian-Lagrangian method for simulation of atomization process.  Comput. Fluids, 176:285–301, 2018.
  * [79] S.A. Orszag.  Analytical theories of turbulence.  J. Fluid Mech., 41:363–386, 1970.
  * [80] P.K. Yeung and S.B. Pope.  Lagrangian statistics from direct numerical simulations of isotropic turbulence.  J. Fluid Mech., 207:531–586, 1989.
  * [81] J.-B. Ritz and J.P. Caltagirone.  A numerical continuous model for the hydrodynamics of fluid particle systems.  Int. J. Numer. Meth. Fluids, 30:1067–1090, 1999.
  * [82] G. Pianet, S. Vincent, J. Leboi, J.P. Caltagirone, and M. Anderhuber.  Simulating compressible gas bubbles with a smooth volume tracking 1-Fluid method.  Int. J. Mult. Flow, 36:273–283, 2010.
  * [83] S. Vincent, J.C. Brandle de Motta, A. Sarthou, J.-L. Estivalezes, O. Simonin, and E. Climent.  A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows.  J. Comput. Phys, 256:582–614, 2014.
  * [84] T. Sayadi, V. Le Chenadec, S. Zaleski, and S. Popinet.  A convergence study of the one-fluid formulation in a phase inversion application at moderate Reynolds and Weber numbers.  to be published in Notes Numer. Fluid Mech. Multidisciplinary Design, 2019.
  * [85] M. Herrmann.  On simulating primary atomization using the refined level set grid method.  Atomization and Sprays, 21(4):283–301, 2011.
  * [86] S. Vincent, M. Tavares, S. Fleau, S. Mimouni, M. Ould-Rouiss, and J.-L. Estivalezes.  A priori filtering and LES modeling of turbulent two-phase flows - application to phase separation.  Comput. Fluids, 176:245–259, 2018.
  * [87] M. Saeedipour, S. Vincent, and S. Pirker.  Large eddy simulation of turbulent interfacial flows using approximate deconvolution model.  Int. J. Multiphase Flow, 112:286–299, 2019.
  * [88] H. Lamb.  Hydrodynamics.  Cambridge University Press, 1932.
  * [89] P. Sagaut and C. Cambon.  Homogeneous turbulence dynamics.  Cambridge University Press, 2008.